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ABSTRACT 


This is a user’s manual for the computer code for partitioning a centralized con- 
troller into decentralized subcontrollers with applicability to Integrated Flight/ 
Propulsion Control (IFPC). Partitioning of a centralized controller into two sub- 
controllers is described and the algorithm on which the code is based is discussed. 
The algorithm uses parameter optimization of a cost function which is described 
here. The major data structures and functions are described. Specific instructions 
are given. The user is led through an example of an IFPC application. 
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1. Introduction 


Large interconnected systems such as the flight /propulsion systems of modern aircraft 
often exhibit significant coupling between the various subsystems. One example of such 
a system is the Short Take-Off and Landing (STOL) aircraft wherein the forces and mo- 
ments generated by the propulsion system provide control and maneuvering capabilities 
for the aircraft at low speeds. This strong coupling suggests that a centralized control 
design be used, however, a centralized controller which is designed for the integrated plant 
considering all the interconnections between the flight and propulsion subsystems may be 
of high order and may be difficult to implement and validate. Specifically, in aircraft design 
it is the responsibility of the engine manufacturer to ensure that the propulsion system 
will provide the desired performance when installed in the aircraft. The engine manu- 
facturer thus needs a separate engine controller to be able to perform extensive testing 
to assure adequate performance and integrity in the presence of operational and safety 
limits. This requirement suggests the need for decentralized implementation of Integrated 
Flight/ Propulsion Control (IFPC) systems. 

One approach to integrated control design which combines aspects of centralized and 
decentralized control design approaches is currently being developed at the NASA Lewis 
Research Center [ 1 ]. This approach consists of first designing a centralized controller, so 
that all subsystem interconnections axe accounted for in the initial design stage, and then 
partitioning the centralized controller into separately implementable decentralized subcon- 
trollers for individual subsystems. Here, partitioning means representing the high-order 
centralized controller with two or more lower order subcontrollers which have input/ output 
intercoupling such that the overall control law obtained on assembling the sub controllers 
closely approximates the input/output behavior of the centralized controller. 

The computer code described in this user’s manual is designed specifically for IFPC 
application and the notation and terminology used here reflects that application. The 
software described here uses a parameter optimization method to match the performance 
of a centralized controller with a partitioned controller consisting of two decentralized 
subcontrollers for the flight and propulsion systems. This matching will be subject to 
certain subsystem design requirements. The structure is shown in Fig. 1.1, with optional 
feedback paths indicated by dotted lines. 

In the decentralized, hierarchical controller partitioning structure shown in Fig. 1.1, 
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Figure 1.1 Controller Partitioning 

the subscripts and superscripts “a” and “e” refer to airframe and propulsion (engine) 
quantities, respectively, and the subscript “c” refers to commanded quantities. The inter- 
face variables z represent propulsion system quantities that affect the airframe, such as 
propulsive forces and moments. The structure is hierarchical in that the airframe (flight) 
controller produces commands for the engine controller via the interface variable (z e<lc ) 
which are tracked by the propulsion subsystem. 

Such a control structure allows the engine manufacturer to evaluate the engine sub- 
system performance independently of the airframe control and to verify that the engine 
subsystem will provide the desired performance when installed in the airframe. In general 
there are practical constraints on the achievable bandwidth of z ea tracking for the engine 
subcontroller. A lower bound on the z ea command tracking bandwidth is based on achiev- 
ing the desired performance for the integrated system, while an upper bound is imposed 
by actuator limits and robustness requirements to high frequency modeling uncertainties. 

The software discussed here refers to the structure described above. The variables are 
named according to the convention given above. The parameters in this optimization pro- 
cess are entries in the state-space representations of the subcontrollers. These parameters 
are bounded so as to maintain subcontroller (open-loop) stability. An assumption made 
in the formulation is that the plant has no direct feedthrough from control inputs, i.e. the 
plant “D” matrix is zero. This simplifies the determination of the cost function and its 
gradient. 

One feature of the software is that the user may separately optimize the airframe 
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controller for a fixed engine controller or optimize the engine controller for a fixed airframe 
controller. The main alternative is to jointly optimize both although separate optimization 
is demonstrated in the example. 

This user’s manual is organized as follows. Section 2 briefly describes the cost function 
which is the objective to be minimized. The partitioning algorithm is described in section 
3. There is also a brief description here of the interrelation among User-defined Functions 
(UDFs) so the user who wishes to change the cost function or partitioning structure will 
know which UDFs must be changed. Section 4 describes the data structures needed for 
carrying out partitioning using the MATRIX x programming language. This section also 
contains a description of the major data structures and variables used. Section 5 has a 
brief outline of the procedure for using the software. Section 6 contains a detailed example 
which exercises the algorithm showing applications of its options. 

Appendix I contain a detailed discussion of the parameterization, the cost function and 
the gradient evaluation as they are implemented in the software. Appendix II contains 
short descriptions of the UDFs which implement the partitioning algorithm. Appendix III 
contains fully-documented source code for partitioning. Appendix IV contains the data 
file, INIT.DAT, and partitioned subcontrollers, [SKA.OPT, SKE.OPT] for the example in 
section 6. 
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2. The Cost Function to be Minimized 


With reference to Figure 1.1, the partitioning problem can be stated as follows: 
Given a centralized controller with transfer matrix K($) and a specifi- 
cation of the partitioning structure of controller inputs and outputs, i.e. 
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where e M = z tac - z«,, so that the closed-loop performance with the sub- 
controllers closely matches that with the centralized controller within the 
requirements of the subsystem. 

The particular subsystem constraint for IFPC application is that the engine subcon- 
troller A' e (s) should have the structure of a command tracking controller for the interface 
variable commands z a0ic . 

The cost function is formulated to reflect the difference between the centralized and 
partitioned controllers. The state space representations of the sub controllers I\ a (s ) and 
K e (s) are parameterized and the cost function is minimized over those parameters denoted 
as a vector p. The formulation of this parameterization is discussed in Appendix I. Sta- 
bilitv robustness may be achieved through the use of optional (user -provided) weighting 
matrices and a normalization function in determining the partitioning cost. Specific details 
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concerning the cost function, the parameters involved and the evaluation of the gradient 
of the cost are contained in Appendix I. 

The cost function /( p) is the sum of the performance cost, /o(p), and an additional 
cost of tracking the airframe-to-engine commands, /i(p), /( p) = /o(p) + /i(p)- 

The performance cost, / 0 , is the i/ 2 norm of the weighted (and possibly normalized) 
difference of transfer matrices for the centralized and partitioned controllers 

OO 

/o(p) = J 
0 


(w o (ju>)(K(ju) - K(p)(ju>))Wi(ju>)^ du ( 2 . 2 ) 


where K is the transfer matrix from the J inputs to the u outputs for the centralized 

controller, and K(s) is the transfer matrix of an “equivalent” centralized controller (having 
the same input /output structure as K) obtained by assembling the partitioned subcon- 
trollers using appropriate plant information. Details of the state space representation for 
K(s) are given in Appendix I. 

We are using the i? 2 norm of the weighted difference between the transfer matrices 
for the centralized controller and the equivalent heir ar chi cally partitioned sub controllers 
as will be described in Appendix I. Since this difference must be strictly proper in order 
to apply this norm, it is reasonable for the D matrices for the centralized and partitioned 
controllers to be the same. Thus it may be desirable to fix the values of D“ a , D* ya , D\ t , 
and D\ yc (as described in Appendix I) to values determined directly by the centralized 
controller. This is one of the options available in “fixing the D-parameters” . 

Wi(juj) and W 0 (juj) are optional input and output weighting matrices, Aperf(w) i s 
an optional scalar normalization function. For example, the weighting 

W i W = G( S )(/ + A'( S )G( S )r I 


has been shown by Dale Enns [ 2 ] to lead to stability robustness for the partitioned system 
provided that the centralized system has this property. Other weighting and normalization 
will be discussed with the example. 

/i (p) is the cost of tracking the z eac command generated by the airframe subcontroller 
for the engine subcontroller. This cost minimizes the difference between the transfer ma- 
trices for the responses to the z a< . command of z ea<; using the partitioned controller and 
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z ea using the centralized controller. 


/i(p) = / £ 
o 1 


A, 

•A^track 


(l|ii..U")-r , '(p)0«)lb) S ‘i" 


(2.3) 


T‘ is the transfer function vector from the airframe commands z Qc to the i lh interface 

cent 

variable with the centralized controller. T* is the i th row of the transfer function 
matrix T from the airframe commands z ttc to the interface variables as commanded by 
the partitioned airframe controller, z e0c , with the partitioned subcontrollers. A, is a scalar 
weighting which determines the influence of / 2 on the total cost and iV TRA CKi(w) are (op- 
tional) scalar normalization functions. || • || 2 denotes the Euclidean norm of the row vector. 
Here, one may use the normalizations Arrack; = I lucent 111 to provide adequate scaling for 
this cost. The parameters A, provide weighting for the contribution of the tracking cost to 
the total cost. It was shown in [ 3 ] that manipulating A, provides an indirect means for 
maintaining reasonable bounds on the z ea c command tracking bandwidth. 

It may be required that the engine subsystem be proper, a condition which would be 
violated if D € eea (described in Appendix I) is nonzero. As a result of the optimization 
process, D \ ta may become large. This possibility is removed by “fixing T>| ea = 0 when 
the option is presented while running the code. 


6 


3. The Partitioning Algorithm 
and its Implementation 

The objective is to minimize the cost /( p) = / 0 (p) + /i(p) as described above where 
the parameters p are certain entries in the state space representation matrices for I\ a (s) 
and K e (s) (denoted SKA and SKE in the code). 

The fixed data used by the algorithm are state-space representations for the plant 
transfer matrix G(s), the centralized controller K($), the (optional) weighting matrices 
W,-($) and W 0 (s) (denoted as SP, SC, SWI and SWO respectively in the code), as well as 
a partitioning structure for the numbers of controller inputs (airframe, MA; and engine, 
ME), numbers of outputs (airframe, KA; and engine, KE), numbers of plant measurements 
(airframe, LA; and engine, LE) and numbers of airframe to engine sub controllers inter- 
face variables (PEA). The control designer may also introduce normalization functions 
(N_PERF and N .TRACK) for the performance and tracking costs. Examples of normal- 
izations are given with the example in Section 6. The user must also enter values of the 
tracking weight parameters A,' which determine the relative contribution of f\ to the total 
cost. 

The algorithm incorporates the Broyden-Fletcher-Shanno-Goldfarb (BFGS) quasi- New- 
ton method to select directions of search for Fletcher’s inaccurate linesearch, see [ 4 ]. This 
iterative method requires the calculation of the combined cost and its gradient for the pa- 
rameters p as referred to above. It uses successive gradients to build up an approximation 
to the inverse Hessian matrix. Moreover, the inaccurate linesearch assures an adequate 
reduction in the cost function at each step without using excessive effort searching for a 
minimum far away from the ultimate solution. In this way, subsequent search steps are 
successively closer to those generated by Newton’s Method and convergence is accelerated 
as the iterations proceed. 

The flow of the parameter optimization algorithm for controller partitioning is shown 
in Figure 3.1. The main steps in the algorithm are: 

1. The initial partitioning is obtained by applying the stepwise procedure described in 
[ 5 ]. Special attention is paid to obtaining reasonably low-order subcontrollers which 
are stable and satisfy the command-tracking requirement. The initial state-space 
representations for the transfer matrices I\ a (s ) and A e (s) are denoted as S_KA and 
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Figure 3.1 Flowchart for Partitioning Optimization Algorithm 

SJKE; note that these axe different from the working representations SKA and SKE 
described next. 

2. The initial partitioning is converted to a ^minimal parameter” form (with state-space 
representations SKA and SKE respectively) and used to generate an initial value of 
the parameter vector P_I. This form is described in Appendix I. 

3. The initial (as well as any subsequent) value of the parameter vector is passed to a 
function which determines the state-space representation for the equivalent partitioned 
controller and calculates the combined cost, / (denoted FP in the code). The gradi- 
ent (denoted as DFDP) is also computed analytically by the procedure described in 
Appendix I. 

4. The BFGS method uses the current gradient in conjunction with previous information 
to generate a direction of search. The Fletcher inaccurate linesearch is carried out using 
the cost and gradient calculated at each parameter vector to predict a new parameter 
vector until one is found which yields a sufficient reduction in both the cost function 
and the size of the gradient. The new point is denoted as PJl. This linesearch is 
constrained so as to maintain stability of the subcontrollers. 
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5. At the end of the lineseaxch, the new parameter vector and cost (PJL1 and FPJ1) are 
compared to the values at the beginning of the linesearch (P J and FP J) as a check on 
convergence. If the maximum change in all the parameters is less than a user-specified 
value and the change in the total cost is less than another value (MAX(P-Il-P-I) < 
EPSILON and ABS(FP Jl-FPJ) < DELTA) then convergence is declared and the 
iteration ceases. If in addition, the maximum absolute value of the partial derivatives 
is less than a user-specified tolerance (ABS(DFDP) < ETA), this is also noted. If the 
number of iterations exceeds ITER or the function value is sufficiently reduced (FP 
< FMIN), then the procedure stops with an appropriate message. If the convergence 
test fails, the algorithm proceeds to update the information used to determine the 
direction of search and to use the most recent cost and gradient values to generate a 
new direction of search and carry out the linesearch via steps 3. and 4. 

6. The output of the algorithm is the state-space representation for sub controllers (de- 
noted SKA_OPT and SKE.OPT) which minimize the cost function /( p) within the 
convergence criteria. 

7. These subcontroller transfer matrices have the same orders n a and n e as the initial 
partitioning. Controller reduction can be performed on these “optimal sub controllers 
and the process of optimization can be repeated on the “new initial partitioning . 

The algorithm is implemented in MATRIXx using a set of functions which are referred 
to in MATRIXx parlance as User-Defined Functions (UDFs). A glossary of variables and 
UDFs follows in Section 4. More complete descriptions of the UDFs appear in Appendix II 
and an annotated version of the code appears in Appendix III. The flow of the MATRIXx 
partitioning code as illustrated in Figure 3.2 follows: 

a. The function START is called with input LAMBDA (required tracking cost weight 
scalar or vector) and STOP (optional stopping criteria vector). 

i. A file (INIT.DAT) is read to acquire the fixed data and initial partitioning (S_KA 
and S_KE) referred to above along with a three dimensional vector (FRQ) which 
gives the left and right end points as well as the number of logarithmically placed 
points in the interval over which numerical integration takes place. 

ii. START calls the routine MODL to put the initial partitioning into an appropriate 
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Figure 3.2 Flowchart for the MATRIX* Partitioning Code 

form to serve as a parameterization for the subcontrollers, i.e. SKA = MODL( S_KA ) 
and SKE = MODL( S _KE ). 

iii. The function LONGCOL transforms these state-space representations for the sub- 
controllers into a parameter vector, i.e. PJ = LONGCOL( SKA, SKE ). This 
function along with its “inverse” MAT ([ SKA, SKE ] = MAT( PJ )) are used 
throughout the code to transform between state-space representations for subcon- 
trollers and a parameter vector. START also determines constants which will be 
used by the other functions and stores them in CONST.DAT . 

iv. START then calls COST to calculate the initial cost (FP) and gradient (DFDP). 

b. START calls PARTITIO with the convergence criteria (STOP) as input. PARTITIO 
is the main routine which 

c. generates the search direction (DI), 

d. calls the function INACCURATE which carries out Fletcher’s inaccurate linesearch 
(bracketing/sectioning) using cost function, F(P J+ALPHA*DI), and gradient values 
generated by the function COST. While INACCURATE is running, the user will see 
displayed first the bracketing interval (AL, ALPR) and function and derivative values 
(FAL, FPAL; FALPR, FPALPR) then the sectioning interval (A, B) and corresponding 


10 


function and derivative values (FA, FPA; FB, FPB). These give the user some sense 
of how rapidly the linesearch is progressing, but they can easily be removed from the 
code if desired. 

e. Convergence/stopping conditions are checked by the function CONVERGE and steps 
c. — e. axe repeated until they are satisfied. During these major iterations, the last 
twenty cost values are plotted. 

f. The output of START is the final partitioning (SKA.OPT, SKE.OPT). At this point, 
the program ends. 

The user may do a posteriori analysis such as order reduction on the subcontroller state- 
space representations (SKA.OPT, SKE.OPT) and start the procedure again with new 
data in the file INIT.DAT. 

If the user wishes to modify the partitioning structure, changes will be necessary in the 
START, COST, LONGCOL and MAT routines. Different constants must be calculated 
and stored by START. Different formulations for the state space representation of K(s ) 
and T(s) must be coded in COST and new formulations for the gradient must be gener- 
ated using the procedure described in Appendix I. Furthermore, the conversions between 
[SKA,SKE] and p by LONGCOL and MAT must be rewritten. 

If a different formulation of the cost function is used, then only the portion of the COST 
UDF where the cost and gradient axe computed must be changed. The algorithm requires 
a gradient with each evaluation of the cost. The new gradient may be the most difficult 
change to make. 
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4. Reference to Major Data Structures and Variables 


Most of the input to the program is provided through the MATRIX x data file INIT.DAT. 
The following data structures are mandatory to the running of the program and must be 
provided in INIT.DAT prior to running the program. 

INIT.DAT — Mandatory Data 

SP, NP — the state space representation and the order for the integrated plant in 
system matrix form 


SC, NSC — the state space system matrix for the centralized controller, and its order. 

S_KA, NS-KA — the state space system for an initial “guess” at the airframe controller, 
and its order. 

S_KE, NSJKE — the state space system for an initial “guess” at the engine controller 
and its order. 

PEA — the number of interface variables from flight controller to the engine controller. 

FRQ — a vector of the form [FRQ(1);FRQ(2);FRQ(3)] where FRQ(3) logarithmically 
placed frequency points over the interval FRQ(l) < u < FRQ(2) are used in the 
numerical integration for determining the costs. The number of points must be 
odd because of the numerical integration rule used. 

The following data structures are optional and may be entered in INIT.DAT if desired. 
Indexing rules of MATRIX x do not allow the index zero or empty vectors. Thus, if 
some of the optioned quantities are absent or have value zero, the code will place dummy 
variables in appropriate matrices and set corresponding size variables to nonzero quantities 
(usually one). 

INIT.DAT — Optional Data — if absent, the indicated default values are set by the code. 

SWI, NWI — state space system for input weighting of the difference between central- 
ized and assembled partitioned controllers, and its order. If absent, the code sets 
SWI = an identity matrix of size (MA -fl ME + LA + LE + 1) and NWI = 1. 
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SWO, NWO — state space system for output weighting of the difference between 
centralized and assembled partitioned controllers, and its order. If absent, the 
code sets SWO = an identity matrix of size (KA + KE + I) and NWO = 1. 

LA — number of integrated plant measurements to the airframe controller. If absent, 
the code sets LA = 1; appropriate zero entries are introduced in SP, SC, and SKA. 
This allows for the case where there are no measurements fed from the plant to the 
airframe controller. 

LE — number of integrated plant measurements to the engine controller. If absent, 
the code sets LE = 1; appropriate zero entries are introduced in SP, SC, and SKE. 
This allows for the case where there are no measurements fed from the plant to the 
engine controller. 

NPERF — normalization vector of size (FRQ(3) x 1) for the performance cost. If 
absent, the code sets NPERF to a vector of ones. 

NTRACK — normalization matrix of size (FRQ(3) x PEA) for the tracking cost. If 
absent, the code sets NTRACK to a matrix of ones. 

STABIL — a necessarily negative parameter which is used to guarantee that all eigen- 
values of the sub controllers have negative real parts for stability. If absent, the 
code sets STABIL = — 10 -9 . 

The major constants used within the code are created by the execution of START and 
stored in the file CONST.DAT. 

MA, ME — number of airframe and engine controller inputs. 

LA, LE — number of integrated plant measurements to airframe and engine controllers. 

KA, KE — number of airframe and engine controller outputs. 

PEA — number of intermediate commands from airframe controller to engine con- 
troller. 

SP, NP — state space system for the integrated plant, with its order. 

S_KA, NS-KA — state space system for the initial airframe controller (J\ a ) in modified 
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modal form (or in original form if the system is to be held fixed) along with its 
order. 

S_KE, NS.KE — state space system for the initial engine controller ( K e ) in modified 
modal form (or in original form if the system is to be held fixed) along with its 

order. 

SK, NK — state space system for centralized controller (A"), with its order. 

SG, NG — state space system for transfer matrix from z ac to (T cen i) using the 
centralized controller, with its order. 

WEIGHT — a vector containing weights to use with Simpson’s integration rule; de- 
pends on FRQ(3) for the number of points at which integration is to take place. 

OMEGA — the vector of logarithmically placed frequency points at which sampling 
is to be done; determined by FRQ. 

AORE — flag to indicate whether the airframe controller (1) or the engine controller 
(2) or neither (0) is held fixed during the optimization process. 

FIXD — flag to indicate whether the ‘D’ matrices are held fixed; if FIXD = 0 none are 
fixed, if FIXD = 1 or 2 the DAA, DAYA, DEE and DEYE are held fixed during the 
optimization process and (if FIXD = 1) then DEEA and DEAE are variables or 
(if FIXD = 2) then DEEA is set to a zero matrix and DEAE is a variable matrix; 
if FIXD = 3 then only DEEA is set to a zero matrix and all the remaining D’s 
are variable matrices. In any case DAA, DAYA, DEE and DE\ E are saved in 
CONST.DAT. 

STABIL — a necessarily negative parameter which is used to guarantee that all eigen- 
values of the sub controllers have negative real parts for stability. If absent, the 
code sets STABIL = -1(T 9 . 

The outputs of interest to the program are kept in the file INTER.DAT, which contains 
a history of the optimization process and information which can be used to restart the 
program (after a crash or after intentionally stopping it) if desired. 
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INTER.DAT contains 


PJ — the last point (parameter vector) that met the inaccurate linesearch minimiza- 
tion criteria, that is, the i th point. P J is a column vector consisting of the successive 
a, /? values in the 2 x 2 blocks of A a and A e , the successive columns of B° a after 
the first, the successive columns of C a , the successive columns of D a correspond- 
ing to the z ea c outputs, the successive columns of B e after the first, the successive 
columns of C| e and the successive columns of D e corresponding to the e ea inputs in 
this order. See Appendix I for a more complete discussion of the parameter vector 
as it is related to the cost function. 

JHO — the complete history of the cost of partitioning through the z lh iteration. 

JH1 — the complete history of the cost of tracking through the z th iteration. 

FH — the complete history of the total cost, FP = FPO + FPl through the i th iteration. 

GRADO — the gradient of the partitioning part of the cost function, /o, at the z th 
iteration. 

GRAD1 — the gradient of the tracking part of the cost function, /o, at the i th iteration. 

GRADI — the total gradient of the cost function at the z lh iteration. 

LAMBDA — the PEAxl vector which weights the contribution of the airframe to 
engine command tracking cost in the toted cost function. 

HI — the inverse Hessian matrix being used during the BFGS optimization process. 

FX — keeps track of which parameters corresponding to A a and A e are at the stability 
bound. 

Refer to the cost section of Appendix I for the structures of SP, SK, SKA, SKE and SG. 
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5. Instructions to the User 


The following instructions are intended to be a quick introduction to the code. They 
do not attempt to explain the details of what is occurring during the execution. More 
thorough documentation of the MATRIXx functions is available in the documented code 
contained in Appendix III. 

• Use MATRIXx to construct the data file INIT.DAT described in the preceding section. 
The initial approximation to the subcontroller state-space matrices S.KA and S_KE 
can be produced by the procedure described in [ 5 ] or may come from a previous 
application of the partitioning software. If any of the optional data are not present in 
INIT.DAT then the code will produce the defaults indicated earlier. 

• Start MATRIXx and type 

DEFINE 1 START. MTX’ 

to activate the code. To execute START you must enter a value of LAMBDA, a weight- 
ing of the tracking cost relative to the total cost. LAMBDA can be a PEA x 1 vector 
whose entries individually weight the z m output responses to the total z a< . inputs. If all 
the weights are to be the same then a scalar may be entered. This variable emphasizes 
the degree to which the tracking cost will affect the total cost. 

Optionally, the vector STOP of stopping conditions can be defined. 

STOP = [EPSL; DELTA; ETA; ITER; FMIN]. 

The program stops if the following criteria are met: 

the maximum change in the parameters MAX(|P J-P Jl|) < EPSL and 


the change in the cost |FP_I-FP-I1| < DELTA and 


the norm of the gradient |DFDP _I| < ETA or 
the number of major iterations I > ITER or 
the cost FP J1 < FMIN 

If STOP is not entered, the following values axe set by the code 
STOP = [EPSL; DELTA; ETA; ITER; FMIN] = [1(T 9 ; 1(T 9 ; 10" 9 ; 100; 0.1]. 
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• Now the program can be executed by entering 

[SKA_OPT,SKE_OPT]=START(LAMBDA, STOP) or 
[SKA_OPT,SKE_OPT]=START(LAMBDA). 

The choice of fixing the engine or airframe system matrices will be presented. If 
neither is to be fixed (the usual choice) enter 0. 

Various options for fixing the D -submatrices are presented. If none are to be fixed 
enter 0 here. 

During the linesearch procedure, the user will see function and derivative values 
which indicate the progress of the search for a reduction in the function value. 
After every major iteration of the program (starting from the second), a MATRIX x 
graph is generated showing the costs (total, partitioning and tracking) for the last 
twenty iterations. 

• When the run ends, the output of START is the final optimized state-space represen- 
tations for the subcontrollers SKA.OPT and SKE.OPT. 

The data file INTER.DAT stores the history of the costs, the most recent parameter 
vector, approximation to the inverse Hessian matrix and gradient. This information 
can be used to see the progress of the algorithm, and to restart the program either 
after a successful termination or after user interruption. 

There is a routine called RESTART which is available for restarting from a system crash 
or an intentional interruption. It uses the data stored in INTER.DAT and CONST.DAT 
and allows the user to define new values of the tracking cost weight A, the stopping criteria 
STOP, and/or a new inverse Hessian approximation. 

• The inputs and outputs of the RESTART routine are similar to those for START. 
There are three alternatives: 

1) enter [SKA_OPT,SKE_OPT]=RESTART(LAMBDANEW) if all that is changed is 
the LAMBDA weighting parameter or 

2) enter [SKA_OPT,SKE_OPT]=RESTART(LAMBDANEW, STOP) if a change is 
made in LAMBDA and/or STOP or 

3) or enter [SKA_OPT,SKE_OPT]=RESTART(LAMBDANEW, STOP, 1) if a restart 
with the identity matrix as the initial approximation to the inverse Hessian is 
desired. Note that you must enter values of LAMBDANEW and STOP even if 
they are the same as the previous LAMBDA and STOP. 
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The restarted program executes in the same way as before, with results stored in 
the file INTER.DAT and the output of the program being the optimized state-space 
representations SKA.OPT and SKE-OPT. 
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6. Example of Controller Partitioning 

STOL Example. 

The controller partitioning software is first applied to a centralized flight /propulsion 
controller for a STOL aircraft as was described in reference [ 6 ]. This controller has the 
form u = iv(s)e with the error vector e consisting of errors, e = [e„, e 9 , em, zeprV ■> * n 
following velocity (u), pitch rate variable (q v = 9 + 0.10), engine fan speed (N 2) and engine 
pressure ratio ( EPR ) commands. The control input vector u consists of rates of change 
of thrust vectoring angle, fuel flow, thrust, reverser port area and nozzle throat area, u = 
[8tv, WF, +78, +8] T . u consists of rates because integrators were appended to the control 
inputs during the process of centralized control design to achieve zero steady-state error for 
step commands. The partitioned airframe and engine controllers are desired to have inputs 
e a = [e„,e g ] T and e e = [e/v 2 i e EP«] T and outputs u a = and u c = [I# ,+78, +8] T 
respectively. The interface variable z m for this example is the single variable FEX, the 
axial thrust generated by the propulsion system. An initial controller partitioning was 
obtained using the procedure discussed in [ 5 ]. 

The numbers of airframe and propulsion subcontroller inputs axe thus (MA=)2 and 
(ME=)2 while the sub controllers have (KA=)1 and (KE=)3 outputs respectively. There is 
(PEA=)1 interface variable and no direct measurements are fed back from the integrated 
plant to the subcontrollers (LA and LE are absent from INIT.DAT since there are no 
measurements). 

State-space matrices for the integrated plant, SP of order (NP=)13, the centralized 
controller, SC of order (NSC=)13 and initial partitioning, S_KA of order (NS_KA=)10 
and S_KE of order (NS_KE=)7 are listed in Appendix IV. 

The optimization is done over the frequency range u> € [0.1, 100] with 41 frequency 
points (FRQ=[0.1; 100; 41]). The frequency weighting 

^,( s ) = G(s)(/ + ^( s )C(s)r 1 

is used to achieve good performance matching as well as stability robustness for the equiva- 
lent controller. The state-space representation for this weighting, SWI of order (NWI=)26, 
is obtained from SP and SC. The tracking normalization, N.TRACK= ||Tc e ntO w )ll 2 i 
used to scale the tracking cost. The state-space representation for Tcent(s) is constructed 
as in expression (1.3) in Appendix I. Since LA and LE are absent, any blocks involving y a 
or y e as either inputs or outputs are omitted. Notice that since PEA=1, Teem is a lx MA 
row vector. Neither output weighting (SWO) nor performance normalization (N-PERF) 
were used. 

All the necessary variables are stored in INIT.DAT . 

The tracking weighting parameter LAMBDA is set to 0.05 and the stopping criteria 
vector is defined as 

STOP=[le-9;le-9;le-9;100;0.l]. 
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After defining the UDF START by entering 

define ’START. MTX’ 

the program is started by entering 

[SKA, SI<E]=START(LAMBDA, STOP). 

Respond to the question concerning fixing the engine or airframe by entering ‘O’ to fix 
neither as shown in Figure 6.1. 

The ‘O’ response to the question concerning the D- matrices as shown in Figure 6.1 
will fix none of these submatrices . The program now begins. 


ENTER 1 TO FIX AIRFRAME, 2 TO FIX ENGINE or 0 FOR NEITHER: 0 

ENTER 1 TO FIX ALL Ds EXCEPT DEAA & DEEA, 2 TO INCLUDE DEEA, 
3 FOR ONLY DEEA, or 0 for NONE: 0 


Figure 6.1 Screen After Responding to Questions 


During a major iteration, the user will see values of AL, FAL, FPAL, ALPR, FALPR, 
A, FA, FPA, B, FB and FPB displayed on the screen. These values result from calculations 
during the lineseaxch as described in Step 4 of Section 3. In particular [AL, ALPR] is the 
interval used in the “bracketing” phase of the linesearch and [A, B] is the interval used 
during “sectioning”. The function values ( and directional derivatives) at the endpoints of 
these intervals are denoted by FAL, FALPR, FA, and FP (respectively FPAL, FPALPR, 
FPA and FPB). The user can follow the progress of the linesearch by viewing the values 
displayed on the screen as in Figure 6.2. 
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AL 


A 

7.5865D-09 


0 . 

FAL 


FA 

4.9981D+01 


4.9737D+01 

FPAL 


FPA 

6.7610D+08 


-6.1169D+08 

ALPR 


B 

0 . 


7.5865D-09 

FALPR 


FB 

4.9737D+01 


4.9981D+01 

FPALPR 


FPB 

-6.1169D+08 


6.7610D+08 

During Bracketing 

During Sectioning 


Figure 6.2 Screen Display During the Linesearch 

After the first linesearch succeeds the user will see a graphical display of the values 
of the total, performance and tracking costs for the previous major iterations (after the 
twentieth, only the last twenty are displayed). A typical screen is displayed in Figure 6.3. 

After 100 iterations, SKA and SKE are returned. The convergence criteria were not all 
met, rather the program stopped because the maximum number of iterations was reached. 
Nonetheless, as will be seen by a posterioiri analysis, the resulting sub controllers exhibited 
good performance and tracking properties. The total cost history is shown in Figure 6.4. 

The performance of the initial controller partitioning is evaluated in comparison with 
that of the centralized controller by comparing closed-loop system response to step com- 
mands in the controlled variables z. 
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Figure 6.3 Screen Display Showing Total, Performance and Tracking Costs 
The responses to q Vc , N2 C and EPR C with the initial as well as the optimized parti 
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Figure 6.4 Cost History for Controller Partitioning Optimization 


tioned subcontrollers were comparable to those with the centralized controller so they are 
not shown here. However, the responses to V c , shown in Fig. 6.5, show considerable 
degradation in terms of increased coupling in the N 2 and EPR responses with the ini- 
tial partitioned sub controllers. This deficiency was overcome by the optimized partitioned 
controllers as can be seen in Fig. 6.5. Note that all the quantities shown in Fig. 6.5 are 
normalized, using scalings discussed in [ 6 ], to allow a direct comparison of the various 
response magnitudes. In addition, the response of FEX (the interface variable) to V c using 
partitioned subcontrollers was also comparable to that using the centralized controller as 
is seen in Fig. 6.6. 

Since the performance with the optimized subcontrollers is found to be acceptable, 
an effort is made to reduce the orders of the subcontrollers. The engine subcontroller is 
reduced to 4 th order by residualization of the three high frequency modes without any 
loss of performance. Through the use of internally balanced reduction techniques [ 7 ], 
the airframe subcontroller is reduced to 6 th order (from the original 10 th order) without 
excessive mismatch in the controller transfer matrix characteristics as is seen from the 


23 


1.5 





Figure 6.5 Closed-Loop System Response to Step Velocity Command for Cen- 
tralized Controller, Initial Partitioning and Optimized Partitioning 
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Figure 6.6 FEX Response to Step Velocity Command for Centralized Controller 
and Optimized Partitioning 

full and reduced order airframe controller singular values comparison in Fig. 6.7. This 
reduced order airframe subcontroller does, however, exhibit deterioration in closed-loop 
performance in the V and q v response comparison plots for a step change in V c as shown 
in Fig. 6.8. 

The reduced order subcontroller state-space matrices are stored as S_KA (order NS_KA=6) 
and SJCE (order NS_KE=4) in INIT.DAT and the program is started by entering 

[SKA, SKE]=START(LAMBDA, STOP). 

One should not use RESTART here since the CONST.DAT file will not contain data related 
to the reduced order subcontrollers. Moreover, since the engine subcontroller is acceptable, 
the optimization should take place over only the airframe subcontroller parameters. When 
requested to enter a value to fix a subcontroller, enter ’2’ to fix the engine subcontroller. 

The program will execute as before, generating optimal SKA and SKE (fixed to the 
initial reduced order S_KE). The response obtained with the optimized reduced order air- 
frame subcontroller for step V c is shown in Fig. 6.8. Note that the optimized subcontroller 
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Figure 6.7 Singular Values of the Airframe Subcontroller for Optimized Parti- 
tioning with A = 0.05 — Full (10) and Reduced (6) Orders 

also provides improved tracking of the velocity command. The state-space matrix for the 
optimized reduced order airframe subcontroller is listed in Appendix IV. 
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Figure 6.8 Closed-Loop System Response to Step Velocity Command for Cen- 
tralized Controller and Reduced Order Partition Subcontrollers — Initial and 
Optimized 


For a discussion of the application of this code to the design of a decentralized controller 
for a Short TakeOff and Vertical Landing (STOVL) aircraft, the user is referred to [ 8 ]. 
This application uses measurements from the integrated plant to the subcontrollers and 
includes more than one interface variable. In addition, several different weighting matrices 
are discussed in [ 8 ]. 
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Appendix I 


Development of the Cost and Its Gradient 


The Parameters. Parameters in the optimization process are certain entries in state-space real- 
izations of K a (s) and K e (s ) as defined in the formula (1.1). The notation M 3 ol is used throughout 
to indicate the matrix M € {A, B,C, D } in the state-space realization of the system transfer matrix 
5 € {c,p,a,e} (c =centralized controller, p =plant, a =airframe subcontroller, e =engine subcon- 
troller) with input i (respectively output o) 6 {p, a,e,ea c ,ea) ( ea c =interface variable commands, 
ea =interface variables). 


K a (s) : 


u Q 

Zen 


C a 
C* a 

ea c a 


D a 

ay<x 

D a 

ea c !/a 


e a 

y Q 


(si -A’)-' 

K'(s): u« = B‘„ B’ t J + (D%„ DU Dl,.)) 

For the purposes of this software description, the corresponding state space matrices are written as 

(I-l) 


®ea 

e e 

ye 



' A a 

B a aa 

1 


SKE = 

’ A e 

B e 

BU 

Bl 1 

SKA = 

C a 

^ aa 

CL a 

_ ea c a 

D a aa 

D ta c a 

D a ayo 

Dea c y a . 

and 

r e 
L We 

^ eea 

D e 

eea 

ee 

DU 

DlyA 


One consideration in choosing a parameterization is to introduce a “minimal’ number of pa- 
rameters in the optimization process. A real canonical form used in [ 9 ] served as the model for 
our parameterization. The subcontroller system dynamics matrices A a and A e are represented as 


0 1 
a 0 


If the order 


block diagonal matrices with two-by-two real companion blocks of the form 
of either A a or A e is odd, there is also one diagonal real entry corresponding to a real eigenvalue. 
In addition, a and 0 are constrained to be negative in order to meet the requirement that sub- 
controllers be stable. It should be noted that this form for the A matrices does not allow for a 
Jordan block structure. However, since the matrices are obtained from a numerical process, it is 
improbable that the “optimal” solution would need such a special structure. 

In addition the first columns of each of the subcontroller input matrices and B e eea are fixed 
at non-zero values determined as follows. State space representations for the initial partitioning 
are not required to be in the canonical form described above. For each of them, a similarity 
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transformation, T, is applied to the initial A matrix so that T AT is in the proper form. If A is 

an n x n matrix then there will be n degrees of freedom in the determination of T . Different T will 
yield the same canonical TAT~ X but different transformed TB and CT' 1 matrices. This implies 
that there are actually n degrees of freedom in the determination of TB and CT' 1 . We select a 
simple T which is nonsingular and compute TB and CT ~ X . We remove the degrees of freedom by 
fixing the n entries in the first column of TB to their values or 10" 9 if the corresponding value is 

zero. 

We are using the H 2 norm of the weighted difference between the transfer matrices for the cen- 
tralized controller and the equivalent heir arch ically partitioned sub controllers as will be described 
below. Since this difference must be strictly proper in order to apply this norm, it is reasonable 
for the D matrices for the centralized and partitioned controllers to be the same. Thus it may 
be desirable to fix the values of D a aa , D a aya , D e ee , and D e eyc to values determined directly by the 
centralized controller. This is one of the options available in “fixing the D-parameters . 


The parameters over which the optimization takes place are then the a and /? entries in the 
block canonical forms, the entries in all but the first columns of the matrices B aa and B eea and all 


the entries in the matrices 


C a 

^afl 

c a 

, u ca c a J 


5 C e %, a and D e eea . The parameter vector will be denoted 
as p € R n where N = n 0 (fc 0 + m a + l a + Pea) + n e (k e + m e + l e + Pea ) + Pea{k a + m e); n , m i k, 
and l refer respectively to the order, number of error inputs, number of outputs, and number of 
direct measurements for a subcontroller and pea refers to the number of interface variables. The 
number of parameters depends not only on the total numbers of controller inputs and outputs and 
interface variables which are fixed but also on the orders of the subcontrollers, n a and n e . There 
is thus a double incentive for keeping these orders low; not only to reduce the complexities of the 
sub controllers but also to accelerate the optimization algorithm whose performance depends on the 
total number of parameters. 


The Cost Function. The cost function is the sum of /o(p) as in (2.2) and /i(p) as in (2.3). 
These involve the transfer matrices K(s), K(s ), T cent (s), and T(s) which are described below. 
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The parameters in the cost function are the entries in SKA and SKE as described in the previous 
section. These parameters along with some fixed transfer matrices are used in the determination 
of the transfer matrices needed in the cost function. 


State-space representations for the centralized controller transfer matrix, K(s), and the plant 
transfer matrix, G(s) = , are given. Those for K(s), T cenl (s) and f(s) (shown in formulas 

(1.2), (1.3) and (1.4) below) are constructed from the state-space representations of K(s), I( a (s), 
K e (s), Gea(s) and <3(s). 


The transfer matrix K(s) which enters into the performance cost term /o(p) depends on K “(s) 
and K e (s), and on the transfer submatrix of the plant from control inputs (u) to interface variables 


(Zee), G ro : z* = C^sl - A p )~ i [B p a B p pe ] 


Ua 

U e 


The block diagram in Fig. 1.1 shows the 


specific interconnections accounted for in this transfer matrix. Note that K (s) has the same inputs 
and outputs as K{s) as described in (2.1). 

A state-space realization for the equivalent partitioned controller K(s) = C(sl - A)~ l B + D 


was shown in [ 10 ] to be S = 


A B 
C D 


where 


A = 

B = 

C = 
D = 


.(B p a cr 


Re r a 

D eea^ ea c 


aa + B P pe Dl 


r a 

eea^ea c 


,{B p a D “ 

Ci a 

D*C1 

. eea ea c 

D a 

D*D* 

eea ea c 


B a aa 

Btea D ta, 


+ B p D e 

aa ' ^pe e 


a 

D a 

eea ea c 


C e 

°ee 


D e ee 


0 

A e 

BUCU 


BU 

BU D 'ee 


-D e eea C p mvJ 
Day* 

Dl'aDU.y. 


-BUaCL, 

(A p - BUD\ ea C p av )\ 


(B p e D‘ 


B a 

ay* 

B'DU „ 

eea ea c y a 

D a + B p D a ) 

eea ea c y* ' pa ay* / 


Dly* 


B\y> 

BV peD\ y * 


( 1 . 2 ) 


The calculation of the tracking cost in (2.3) requires two transfer matrices, T cen t(s) and T(s). 
The norms of the rows of their difference measures the differences in response of the various possible 
interface variables to airframe commands. Figure 1.2 is the block diagram for T cen t(s) = C cent (s7 — 
^cent^-ijgcent _j_ £>cent ^ t h e c i ose( j i DO p transfer matrix from airframe commands, z ac , to interface 

variable quantities, z^, produced by the engine using the centralized controller. 
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^pea 
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JTp” 

l u pye 


-Ye 


Figure 1.1 Partitioned Controller 
Note that in Figures 1.1, 1.2 and 1.3 a block of the form 


A a 


represents the block 



As is easily seen from Figure 1.2 the state-space representation for the centralized controller 
transfer matrix T ce n t(s) be written in terms of the submatrices in the state-space representations 
for K(s) and G(s ), 


A cent = 


^gcent _ 


(B; a c< 


ac + BpeCZc) 


(-B' ca CZ P - Bl t C% + B^Cl P + B‘ cy Cl P 

~ + Dt'CL + D c 


B c 

^ ca 

B v D c 

. ‘ Ly pa J - y aa. . 


A^ + B; a (-Dt a CZ p + D c a 
+ BlX~D\ K Cl p + D 

C cenl = [0 C£ p ] 


rv \ 

ayc^ycpJ 


rv 

QVo y«P 1 ^ ae'-'ep 

i,.cip + DC "C r , r + d;.. cr. 


ey c ~y*v' 


'jcent 


= o. 


(1.3) 


Figure 1.3 is the block diagram for T(s) = C(sl - A)~ l B + D , the closed loop transfer matrix 
from airframe commands, z ac , to interface variable quantities produced by the airframe subcon- 
troller, z ea c using the partitioned controller. 
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Figure 1.2 z m response to z ac using the centralized controller 
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Figure 1.3 z ea<: response to z ac using the partitioned controller 


33 














The state space matrices for T(s) = C (si — A) 1 B + D can be written in terms of the st ate- space 
representations foT K '“(s), A’ e (s) and G(s) as 


T = 

A B 
C D_ 

where 





A a 

0 

-B a C v + B a C p 1 

^aa^ap ' ctVa VaP 

A = 


B' eea C^ a 

A e 

A e p 


ftp r “ -I- B v D e C a 

\- n pa'-'aa ' ^pe^eea^ ea c a 

B? e C' ee 

A pp 

B = 

B a aa 

Btea^a 

B p pe D a aa + Bt'D^D^ J 



c = [c^ ea 

0 -Di..c; r + 

] 


D = (D a mca ) 


(1.4) 


where A„ = (- B’,C! P - C' tn Cl, p - + B’„D^ t ,SL e + “ d *rr = ( A ' ~ 

Bt.DUC{, - - B’„D im Q, - B;,D'«C! r + C% + CJ.„ + 

b;,duci t ). 

It may be required that the engine subsystem be strictly proper, a condition which would be 
violated if D e eea is nonzero. As a result of the optimization process, D e eea may become large. This 
possibility is removed by “fixing D e eea = 0” when the option is presented while running the code. 


The total cost is evaluated for a particular parameter vector p (corresponding to particular 
SKA and SKE) by applying Simpson’s Rule for numerical integration to / 0 (p) + /i(p) over a 
user defined logarithmically spaced frequency interval [wi,wa]. The expressions given above for 
the state-space representations of K(s),R(s),T cen t(s) and T(s) are used for calculating /o(p) &nd 
/j(p) according to formulas (2.2) and (2.3). 


The Gradient of the Cost. 

The performance cost /o(p) was defined in (2.2) as 


/o(p) = / 



A'(p))Wi)* (W 9 (K 


K(p)Wi) 


dw 
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where the explicit dependence of the integrand on u> is suppressed for convenience. 

Only K( p) depends on p so the derivative of /o(p) with respect to a parameter p of p is 


fl/o(p) 

dp 


= / ( Wo{K ~ K{p )Wi Y (' Wo{K ~ K ^ Wi ) 

0 

+ fw o (K - /r£p))Wi)" Y (w o (K - A>))Wi)]du; 


= - 2 / n^ u - K(p))W -y w ° dJ TT w ' 


du; 


Since Jif(p)(s) = C(sl — A ) 1 B + the product rule for differentiation implies that 

d J^M 2)-. b + c(»; - 2)- i( 5 / - 2)- * + c ( ./ - 2)-» f + f 


dp 


dp 


dp 
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= [C(sI-A )-' 7]^ 


'a b' 

'(. sI-A)~ l B ' 

C D. 
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Thus, 

d/o(p) 

dp 


OO 

= - 2K Vw r [(^ ( *~^ ))w 4 


w.[cow-2)-' /]^ 


4 5 

C 2? 


(ju;7 - 


/) ’ 5 ] IV,] 


du> 


Using abbreviations for the system matrix given in (1.2) well as for the terms on the left and 
right sides of the partial, 


5 = 


4 & 

C D 


L{ w) = — — (w / 0 (7v - A'(p))Wi) V. [C(ju>I - I)" 1 7] 


R(u) = 


•^PERF 

(jul- A )- 1 B 

I 


Wi 


allows the derivative to be written as 

df o(p) 


dp 


= -2 Re J tr X(u>)^-7£(ci>) 


o 

2i-th 


du>. 


Any particular parameter p of p is some jk lh entry of some submatrix, denoted M 0 \, relating 
an input i to an output o of SKA (or SKE) as described in (1.1). Furthermore, p (as well as the 


entire submatrix Ad 0 i) occurs in one or more blocks of 5 = 


4 B 
C D 


as was described in (1.2). 
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S can be thought of as consisting of blocks aligned in “columns” corresponding to the “inputs” 


x 0 , x e , x p , e Q , e e , y a , y e and in “rows” corresponding to the “outputs” x a , x e , x p , u a , u e . Each 
block of S is denoted as Bq i where I is one of the “inputs” above and O is one of the “outputs”. 
Let Eo denote a column block matrix with the same number of rows as S , the same number of 
columns as the dimension of the “output” O, with an identity matrix in the rows corresponding 
to “output” O, and with zeros elsewhere. Let E\ T be a similar row block matrix corresponding to 
the “input” I, 


E o = 


Ej T = [0 • • • 0 I 0 


0 ], 


Using this notation, S can be written as 


S = EoBoiEi T . 

o i 

In the partial derivative of S with respect to the p = M 0 \ jk , every block of S is zero except 
for the blocks containing p- Denote such a block as Boi = Eoo-M 0 iTZ-\i-> where Coo denotes the 
factor to the left of (if one exists, otherwise an identity) and Tl\i denotes the factor to the 
right ( similarly, an identity if M 0 \ is the rightmost factor); define both £ 0o = 0 and TZ- a = 0 if 
Boi does not contain M Q \ ■ The partial of S with respect to p thus contains a term of the form 
Coo e j e k' r 'R-ii in the same place as the block Boi- These contributions can be written as 

^ = EoCo'tj'JKilES 


Thus the partial of the partitioning cost can be expressed as 


d/o(p) 

dp 



EoBoo^i^k T E,\iEi T R{v) 


du. 
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Note that L(u)E 0 C 0o ej is a column vector and ejlln Ei T R{u) is a row vector. Furthermore, 
the trace of a column times a row is the dot product of the row and the column. Therefore, 

OO 

= -2 R e/^E W^iEi T R{^)mu)EoCoo^)} dw. 


o i 


To generate the partials with respect to all the parameters in M ol simultaneously, let j and k 
vary over the row and column indicies respectively of the submatrix. Notice that varying the row 
index j selects the j th column of the product L(u)EoC 0o , whereas varying the column index k 
selects the k lh row of HiiEi T R(u>). 

A matrix containing the partial derivatives with respect to the entries of M Q \ located in each 
entry’s proper spot is thus obtained by replacing ej and e k by identity matrices and transposing 
the result. 

OO 

= -2 R e/EE [(KnE J T R(u>))(L(u)E 0 Eoo)] T du, 

0 o 1 

Note that only the terms in R(u)L(u>) depend on u and thus the integration can be rewritten to 
yield 

fl/o(p) . 

dM 0 \ 

^ * L( 

Finally, build two matrices denoted an< ^ dSKE ° ^ s ^ a ^ es an ^ SKE (de- 

scribed in (1.1)) respectively containing the partial derivatives of /o(p) with respect to the parame- 
ters in SKA (respectively SKE) in the same positions as the corresponding parameters would occur. 
This is done so that in the software the gradient vector can be produced from an( * 

by a call of the function LONGCOL (the same function which produces p from SKA and SKE). 

To build these matrices, define the “row” and “column” block matrices E a \ T and E a 0 relative 
to the “inputs” i € {x a ,e a ,y Q } and “outputs” o € {x e ,u a ,Zea c } for the airframe controller state- 
space system matrix SKA. As before E a Q is a block column matrix with as many rows as SKA, 
with an identity matrix of size equal to the dimension of the output o in the rows corresponding 
to o and with zeros in the remaining rows. The matrices column block matrix E e 0 and the row 
block matrices E a i T , and E e \ T are similarly defined. If we again denote the submatrices of SKA 


2 Re 

E E 

OO 

I [R(u)L{u)]du 
J 

EoCoo 


O I 

Lo J 



37 


as M a oi and those of SKE as M e Q i, we can write 


SKA = E E i T Md SKE = E E E'vM'eiE-?. 

o i 0 » 

Notice that, as with 5, the column and row block matrices merely position each M 0 » properly. 
Thus, replacing M OI by the corresponding block in these formulas gives the desired results 

d/o(p) _ pa d/o(P) pa.T 

dSKA'^4- °dM‘ 0 i * 

O 1 


and 


^/o(p)_v~'Y^ pe dfojp} pe t 
dSKE °dM e ■ ‘ 


dfo(p) . 

Using the expression from above for Q ■ gives 


dM a , 


df o(p) 

dSKA 


= - 2Re EE E ’ 


EE 


o i 


00 

1 [R{u})L(u>)}duj 
J 

JSo^Oo 

0 J 

_ 


= —2 Re 


EE EE E °' K “ £ > : 


o i O I 
T 


oo 


1 [R(ui)T(w)]dui 

E 0 CooE a 0 T 

J 

Lo J 

- 


£V 

T 


The terms of the form E* x HnE? are independent of “outputs” o and O whereas those of the form 
EoEo 0 E a 0 T a re independent of “inputs” i and I. Therefore the sums can be rearranged 


d/o(p) _ 


dSKA 


= -2 Re 


(^2^2E a iTZiiEi T ^j ^J[R(u)L{ w)]dwj ^££o£°o£ a o T 


The sum £ E represents a matrix with “inputs” I, the “inputs” of 5, and with “out- 

i I 

puts” i, the “inputs” for SKA. This matrix has the submatrix 72.ii as the block in the rows corre- 
sponding to i and the columns corresponding to I; denote this matrix as R“. Similarly the matrix 

L“ = £ Y, EoEo 0 E a 0 T contains the submatrix £oo in its block with rows corresponding to the 
o o 

5 “output” O and with columns corresponding to SKA “outputs” o. 

, . d/ 0 (p) 

The following simple procedure can thus be used to determine 
1. For each block Z?oi of 5 containing a submatrix M a 0 \ from SKA, determine the left and 
right factors (or identity) £oo and 72-ii of -M a 0 i- 
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2. Enter £ 0o in the appropriate block of L° and TZn in the appropriate block of R“ 

OO 

3. Calculate J R{u)L(uj) du> where 


L(u) = ^4^ (W 0 (K - K(p))W,yW 0 [C(ju>I - A r 1 I] 

(ju)J - A)~ l B 

I 


R(u) = 


Wi 


4. Form 


dfoil P) 
dSKA 


= -2 Re 


R° 


r oo 

J [R(u)L(u>)] dc 


L° 


2 Re 

R£ 1 

J[R(u)L(u)}du] 

1 L ' 



\o / 



By a similar procedure 

d/o(p) . 

3SKE ' 

where L e and R e contain the left and right factors respectively of M e 0 i terms appearing in 5. 
It is easy to use the representation of 5" as given in (1.2) to calculate 


R“ = 









■I 

0 

0 - 

'I 

0 

0 

0 

0 

0 

0‘ 


0 

0 

B e 

^ eea 

0 

0 

0 

I 

0 

0 

0 

and L° = 

0 

B p 

^pa 

B p D € 

eea 

0 

0 

0 

0 

0 

I 

0 


0 

I 

0 








.0 

0 

D e 

eea J 


and 






0 

0 

0 

D a 

o o 
1 


-o 

0 - 

- / 
c e 

'-'eea 

0 

0 

0 

-c* 

eap 

0 

D a 

ea c a 

and L e = 

I 

0 

0 

B* 

0 

0 

0 

0 

I 

0 

0 

1 o o 

J 

pe 

0 

/ . 

. 0 

0 

0 

0 

0 

0 

I. 



R e = 


The Tracking Cost 
The tracking cost was given in (2.3) as 


du> 


/,(p) ” /£ aWkM r ’ (p)lb -* 

0 1 

where the state space representations for T cent (s) and T(s ) are given in (1.3) and (1.4) respectively. 
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This function can be put into a form which is similar to /o(p) so the same procedure for calcu- 
lating the gradient applies. Each row of the difference is normalized before the sum of squares is cal- 
culated. Note that the multiplication by non-negative A, and division by non-negative N 'track ( w ) 
could also be absorbed into the normalized sum of squares by the use of their square roots. Multi- 
plication (or division) of rows by factors can also be achieved by multiplication from the left by a 
diagonal matrix. The sum of squares of the row norms in the resulting product is the same as the 
sum of squares of all the entries or the i7 2 norm of the diagonal weighted difference. In this case 



As before, denote 

L t (u) = 


diag 


A< 


(Tent - f (p)) 


-^TRACK , , 

J nuLm) [5(juI ~ S) " 


TRACK 


(«) 


Rt{u) = 


(jul- A)- 1 B 

I 


and apply the same procedure as earlier to write the partial derivatives as 


d/i(p) 

3SKA 


= -2 Re 


7 0 0 0 

0 0 -C% I 

0 0 Cl, 0 


oo 

J R t (u)Lt(u>)du> 


r/ 0 
0 0 
0 B v B v D e 

u ■ u pa ^pe-^eea 

oo i 


B* 

eea 


and 

d/i(p) 

dSKE 




■ 0 

7 

0 

0 

OO 

= -2 Re 


c e 

* -'em 
0 

0 

0 

~rv _ n e C p + D a C v 
-C v 

ep 

0 D a mca 
0 

J R t (u>)Lt(u>)du> 



. 0 

0 

c p 

0 

0 


I 0 
0 B\ 
L0 0 


This completes the discussion of the cost function and its gradient as implemented in this 
software. If the user wishes to compute cost functions involving the H 2 norm as used here, the user 
must apply the procedure described above to determine the gradient. 
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Appendix II 


Short descriptions of User-Defined Functions 


START in START.MTX 

[SKA_OPT,SKA_OPT] = START(LAMBDA, STOP) 

Initializes information from data files, puts subcontrollers into modified modal form, constructs or 
initializes some matrices required for the evaluation of the cost function. Calls COST to initialize 
the costs and gradients and then calls PARTITION to perform the optimization. 

Input 

LAMBDA — weighting for the contribution of the tracking cost to the total cost. 

STOP (optional) — vector of stopping conditions. 

INIT.DAT — data file containing initial information. 

Output 

SKA.OPT — the state-space representation of the optimized airframe subcontroller. 
SKE.OPT — the state-space representation of the optimized engine sub controller. 
CONST.DAT — data file containing constants used by other UDFs. 

PAR.DAT — used to store the stopping conditions in STOP. 


RESTART in RESTART.MTX 

[SKA_OPT,SKA_OPT] = RESTART(LAMBDANEW, STOP, NEWH) 

Restarts the program using the data available in INTER.DAT and CONST.DAT. Calls COST to 
initialize the costs and gradients, and then calls PARTITION to optimize. 

Input 

LAMBDANEW — a new value of LAMBDA may allow a different emphasis on the tracking 
cost relative to the total cost. 
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STOP (optional) — the vector of stopping conditions, same as in the START routine. It may 
be redefined, perhaps to allow more stringent conditions or more iterations. 

NEWH (optional) — a flag whose presence indicates the desire for restart with identity inverse 
Hessian. If restart with the current Hessian matrix is desired, no value should be passed. 
CONST.DAT — data file containing constants pertaining to the plant and global controller 
being partitioned. 

INTER.DAT — the data file in which intermediate results from previous iterations are stored. 
Output 

SKA.OPT — the state-space representation of the optimized airframe sub controller. 
SKE.OPT — the state-space representation of the optimized engine subcontroller. 

PAR.DAT — the data file in which the stopping conditions STOP are saved. 

INTER.DAT — the data file in which intermediate results are stored. 

PARTITION in PARTITIO.MTX 
PJ = PARTITION(STOP) 

This is the main routine and does the optimization. It iterates till some convergence or stopping 
conditions are met. After each iteration a graph is plotted showing the change in costs. PARTI- 
TION implements the Broyden-Fletcher-Goldfarb-Shanno method of determining a search direction 
for minimization. It then calls the function INACCURATE to implement Fletcher’s inaccurate line 
search in that direction. The function CONVERGE is called to see if convergence or stopping 
conditions are met; if not, the iteration in PARTITION is repeated. 

Input 

STOP — the vector of stopping conditions, used to check convergence. 

CONST.DAT — the data file containing the constants of the program. 

INTER.DAT — the data file containing the intermediate results. 

Output 
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P J — the vector of parameters, after the i th optimization step. 

INTER.DAT — the data file containing the intermediate results, where the results of the 
optimization are stored. 


INACCURATE in LINESRCH.MTX 

[ALI, FAL, FPO, FP1, DFDPAL, DFPO, DFPl, FLAG] 

= INACCURATE (X, FZ, DFDPZ, D, FMIN, DELTAF, LAMBDA, ALPHAMAX) 

This function performs Fletcher’s inaccurate line search as part of the unconstrained optimization 
performed by the function PARTITION. It calls the function COST to get the cost and gradient for 
the cost function evaluated at the i lh set of parameters. Effectively, INACCURATE seeks a point 
where a sufficient decrease in both the function value and directional derivative have occurred. 

Input 

X — the current point (parameter vector), before the linesearch. 

FZ — the total cost function evaluated at X. 

DFDPZ — the gradient of the cost function evaluated at X. 

D — the direction vector which the PARTITION function has chosen to perform the line search. 
FMIN — the minimum value for the cost function. If the cost falls below this value, the 
function will terminate. 

DELTAF — the estimated change in cost. 

LAMBDA — the weighting for the tracking cost. 

ALPHAMAX — maximum alpha, set by stability constraints. 

Output 

ALI — the alpha value that yields sufficient reduction in F(X + alpha x D) for the i th iteration. 
FAL — the total cost function evaluated at X + AL x D. 

FPO — the partitioning cost evaluated at X + AL x D. 
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FP1 — the tracking cost evaluated at X + AL x D. 

DFDPAL — the gradient of the total cost function evaluated at X + AL x D. 

DFPO — the gradient of the partitioning cost at X + AL x D. 

DFP1 — the gradient of the tracking cost at X + AL x D. 

FLAG an output flag signalling the condition with which INACCURATE completed. 

0 — solution found, no problems. 

1 — solution found, objective value less than FMIN. 

2 — new point found, but backeting step converged to right endpoint without satisfying 

Armijo/Goldstein conditions. 


COST in COST.MTX 

[DF, DFO, DF1, FP, FPO, FP1] = COST (P, LAMBDA) 

Evaluates the cost and gradient at the current point using the cost function described in Appendix I. 
The state space representations for K(s) and f (s) are constructed based on the current value of the 
parameter P. The contributions of these to the costs /o and f\ aTe calculated at each frequency point 
in OMEGA and are added according to Simpson’s Rule for numerical integration. Furthermore, the 
contributions to the gradient as described in Appendix I aTe also accumulated for each frequency 
in OMEGA. 

Input 

P — the current point (parameter vector). 

LAMBDA — the weighting for the tracking cost. 

CONST.DAT — the data file containing the constants of the program. 

Output 

DF, DFO, DF1 — the gradients of the total cost function, the partitioning cost function and 
the tracking cost function respectively. 

FP, FPO, FP1 — the values of the total cost function, the partitioning cost function and the 
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tracking cost function respectively. 


CONVERGE in CONVERGE.MTX 

ANSWER = CONVERGE (XJ1, FJ1, XJ, FJ, GRAD, EPS, DEL, ETA) 

This function checks if the convergence condition are met and returns a flag ANSWER defined as 

{ 0 if MAX |XJ1-XJ| >EPS or |FJ1-FJ| >DEL 

1 if MAX |XJ1-XJ| <EPS and |FJ1-FJ| <DEL but MAX |GRAD| >ETA 
2 if MAX |X_L1-XJ| <EPS, |FJ1-FJ| <DEL and MAX |GRAD| <ETA 

Input 

XJ1 — the value of the variable (parameter vector) at the (t + l) th iteration. 

F_I1 — the value of the cost function at the (i + l) th iteration. 

XJ — the value of the variable (parameter vector) at the t th iteration. 

F J — the value of the cost function at the i th iteration. 

GRAD — the gradient of the cost function at the (i + l) th iteration. 

EPS — the absolute tolerance for X differences. 

DEL — the absolute tolerence for cost function differences. 

ETA — the absolute tolerence for the gradient. 

Output 

ANSWER — a flag indicating the state of convergence. 

0 — no convergence. 

1 — variable values converge, function values converge. 

2 — variable values, function values, and gradient converge. 


MODL in MODL.MTX 
[SM]=MODL(S,NS) 

This function takes a system matrix and puts it in to the form where A has 2x2 companion matrix 
blocks whose first rows are [0 1] and whose second rows are [a bj. The transformation matrix is 
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normalized by requiring that all nonzero entries in the first column of the B matrix remain fixed 
(zero entries are set at 10 -9 ). 

Input 

S — the system matrix. 

NS — the order of the input matrix S. 

Output 

SM — the system matrix in modified modal form, still having order NS. 


MAT in PARMAT.MTX 
[SKA,SKE] = MAT(P) 

Creates the system matrices foT the airframe and the engine controllers from the parameter vector 
(the variable over which the optimization process is being performed). In case one of the subsystems 
is fixed, the corresponding initial system matrix is loaded from CONST.DAT and returned as SKA 

or SKE. 

Input 

P — the parameter vector. 

CONST.DAT — the data file of constants. 

Output 

SKA — the system matrix for the airframe controller. 

SKE — the system matrix for the engine controller. 


LONGCOL in PARVEC.MTX 
[P] = LONGCOL(SKA,SKE) 

This function generates the parameter vector from the system matrices for the airframe and the 
engine. In case one of the sub controllers is fixed, the parameter vector corresponds only to the 
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parameters in the other subcontroller. 


Input 

SKA — the airframe controller. 

SKE — the engine controller. 

CONST.DAT — the data file of constants. 

Output 

P — the long column vector of the parameters. 


Z in ZERO.MTX 
[ZER] = Z(NROW, NCOL) 

Constructs a matrix of zeros of size NROW x NCOL. 

Input 

NROW, NCOL — the row and column size of the desired zero matrix. 
Output 

ZER — the generated zero matrix. 


47 


Appendix III 

Controller Partitioning Code 
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convcrge.mtx 

// answer-converge (x_il, f _ll,x_i, f_i, grd, epsilon, delta, eta, fx) 

II 

// The function CONVERGE located in the file CONVERGE. MTX is called by 
// PARTITION to test the convergence of the optimization algorithm according 
// to the following criteria: 

// Condition X: | xil(j)-xl(j) | < epsllon(j) for j-l..n where dim(x)-n x 1 

// whether the parameter vectors are converging 

// Condition f: | f(xil)-f(xi) | < delta 

// whether the function values (costs) are converging 

// Condition G: |grd(J)| < eta for j-l..n where dlm(x)-n x 1 

// whether the gradient Is approaching some low value (approx 0) 

// The result of the convergence check are returned in the flag answer. 

// The codes used for the answer are: 

// answer - 1 if X and f are true and G is false. 

// answer - 2 if X,f and G are true 

// answer - 0 if either X or f is false 

// Input: 

// x_ll - (real,vec) value of variable at iteration (i + 1) . 

// x_i - (real,vec) value of variable at Iteration (i). 

// f_i and f_il - (real, seal) values of function at successive iterations 

// g_i and g_il - (rcal,vec) gradients at successive iterations 

// epsilon - (real,vec) vector of the absolute tolerances defined 
// for each of the elements of the x 

// delta - [real, seal) absolute tolerance for the objective values. 

// eta - (real, seal) absolute tolerance for gradient 

// fx - [real,vecl tracks which of the A matrices is at its stability bound 
// Output: 

// answer - (int,scal) flag that will show whether or not the 
// optimization procedure has converged according to the 

// user-defined criteria. 

// 

answer-0; 

(row, col) -size (grd) ; 

if abs (f_il-f_l) > delta then retf; 
d i f f-abs (x_il-x_i ) ; 
for i-l:row; ... 

if dlff (1)> epsilon (1) , retf;... 
end; 

answer-1; 

for l-l:row; ... 

if f x (1) -0, if grd (1 ) >eta, retf; ... 
end;end; 

answer-2; retf; 


retf 

// Created: 01/24/90 

// Programmer: Steven Ims 

// Revised by Phil Schmidt and students Nader Kamrani and Brian Holawecky at 
// The University of Akron with the support of NASA Lewis Research Center 
// under grant NAG-3-1146. 
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........ ................. 

li^MH 

I0:17:3X j 

// (Df , DfO, Df 1 , fp, fpO, fpl) -cost (p, lambda) 

// 

// The function COST in file COST.MTX 

// computes the cost and gradient of the cost function at the current 
// point (parameter vector). The partitioning and tracking costs and 
// their gradients are returned seperately. 

// Input: 

// p - the parameter vector 

// lambda - the weighting applied to the tracking cost In the total 
// cost function 

// const.dat - datafile of constants 
// Output: 

// Df - gradient of the total cost function 

// DfO - gradient of the partitioning cost function 

// Df 1 - gradient of the tracking cost function 

// fp - the total cost at point p (current parameter vector) 

// fpO - value of the partitioning cost function 

// fpl - value of the tracking cost function 

// 

load ' const .dat ' la le ka ma ke me pea nka nke sp sk np nk frq omega ... 
sg ng swl nwi swo nwo weight stabll aore Nperf Ntrack; 


//Generate ska and ske controllers from the parameter vector p 
( ska, ske) -mat (p) ; 


S\ 

D 


//Note some of 

mt-ma +me; 

lt-la+le; 

mtot-mt+lt; 

kt-ka+ke; 

nt»np+nka+nke; 


the sizes of matrices, numbers of Inputs and outputs 
//ma,me - airframe, engine Inputs 

//la,le - numbers of airframe and engine y feedbacks 

//mtot - total numbers of inputs to the centralllzed controller 

//ka, ke - airframe and engine outputs 

//np, nka, nke - orders of plant , a 1 rf rame, eng ine controllers 


// Split the plant controller sp and find the sizes of each piece 
// using the stored data In lnlt.dat. Get the pieces from each submatrix M 
// which connect input 1 to output o, that Is Mol; Bpa, Bpe, Cap, Cyep 
// Some of the Indices Into the rows and columns used repeatedly are 
// precalculated to save processing, le, mapl-ma+1 etc 


(ap,bp f cp,Df j -split (sp, np) ; 

(mp ntmp) -size (cp) ; 

mapl-ma+1 ;mame-ma+me; 

mpea-ma+me+pea; 

kapl-ka+1; kp-ka+ke; 

bpa-bp(: , ( 1 : k a ) ) ; 

bpe -bp ( : , (kapl :kp) ) ; 

cap-cp ( ( 1 :ma ) , : ) ; 

ceap*cp ( (mapl :ma+pea ) , : ) ; 

cep-cp ( ( (ma+pea+1) :mpea) , :) ; 

cyap-cp(( (mpea+1) : (mpea+la) ), 

cyep-cp(( (mpea+la+1) : (mp) ), :) , 


//split sp 

//the ap: xp -> xp and does not need 
// to be further broken up 
//the dp submatrix of the plant Is 
//kept zero 

//part of bp: ua -> xp 
//part of bp: ue -> xp 
//part of cp: xp -> za 
//part of cp: xp -> zeap 
//part of cp: xp -> ze 
) ; //part of cp: xp -> ya 
//part of cp: xp -> ye 


// From the airframe controller ska, get Baa, Baya , Caa , Ceaa, 

// Daa, Daya , Deaa, Deaya In similar manner to the way that the plant 
// matrix was split up 

(aa , ba, ca, da ) -spl It (ska , nka) ; //split the airframe controller 
(katot ntmp) -si ze (ca) ; //part of 


baa-ba ( : , (1 :ma) ) ; 

// 

ba: 

ea -> 

xa 

baya-ba ( : , ( (ma + 1) : (ma + la) ) ) ; 

// 

ba : 

ya -> 

xa 

caa-ca ( ( 1 :ka ) , : ) ; 

// 

ca : 

xa -> 

ua 

ceaa-ca (( kapl : katot ) , : ) ; 

// 

ca : 

xa -> 

zea 

daa -da ( ( 1 :ka) , (l:ma j) ; 

// 

da : 

ea -> 

ua 

daya -da ((l:ka) , (mapl: (ma + la) ) ) , 

; n 

da : 

ya -> 

ua 

deaa -da ( l kapl : katot ) , (1 :ma) ) ; 

// 

da: 

ea -> 

zea 


deaya-da ( (kapl : katot ), (mapl : (ma + la) )) ; // da: ya -> zea 


// From the engine subcontroller, get Beea, Bee, Beye,Cee, Deeam, Deem, Deye 
(ae, be, ce, del-split (ske, nke); //split the engine controller 
(ntmp metot 1 -si ze (be) ; //part of 


beea -be ( : , ( 1 : pea ) ) ; 

// 

be: 

zea -> xe 

bee-be ( : , ( (pea+1) : (me+pea) ) ) ; 

// 

be: 

ee -> xe 

beye-be (: , ( (me+pea + 1) :metot ) ) ; 

// 

be: 

ye -> xe 

cee-ce; 

// 

ce: 

xe -> ue 

deea-de ( : , ( 1 : pea ) ) ; 

// 

de: 

zea -> ue 

dee-de ( : , ( (pea+1) : (me + pea) 1 ) ; 

// 

de: 

ee -> ue 

deye-de (:, ( (me+pea+1) :metot)) ; 

// 

de: 

ye -> ue 

//PARTITIONING COST MATRICES 

// Get the Centralized controller K 

for 

the performance cost as 

// Input data originally read 

from 

lnlt 

.dat 


(ak,bk, ck,dk) - split (sk , nk) ; 

// Form the Equivalent System tilde (K) using Partitioned Subcontrollers 
// from inputs: Ee,Ea,Y to outputs: Ua,Ue. 
at 0 -(aa, 0 *ones(nka,nke) , 0 ‘ones (nka, np) ; .. 
beea* ceaa, ae, -beea*ceap; . . 

bpa*caa+bpe*deea*ceaa ,bpe*ce, (ap-bpe*deea *ceap) ); 
btO- (baa, 0* ones (nka, me) , baya, 0* ones (nka, le) ; . . 
beea* deaa, bee, beea* deaya, beye; . . 

bpa * daa +bpe* deea* deaa, bpe* dee, bpe* dee a * deaya + bpa* daya , bpe* deye) ; 
ct0-(caa, 0*ones (ka, nke) ,0*ones (ka, np) ;deea*ceaa,ce, (-deea * ceap) ) , 
dtO- ( daa, 0* ones (ka, me) , daya, 0*ones (ka, le) ; . . 
deea* deaa, dee, deea* deaya, deye) ; 

//TRACKING COST MATRICES 

// Form Tcent from input: Zac to output: Zeac for 

// coirjnand tracking using Centralized Controller. This was created in 
// start. mtx because it only needs to be defined once (does not vary 
// with the parameter) 

(ag,bg, cg,dg) -spl It (sg, ng) ; 


// Generate the state-space representation tllde(T) from input:Zac 
// to output: Zeac, command tracking using Partitioned Controller. 
atc-(aa,0*ones(nka,nke) , (-baa*cap+baya*cyap) ; . . . 
beea *ceaa, ae, (-l*bee*cep-beea*ceap-beea*deaa*cap ... 
+beea*deaya*cyap+beye*cyep) ; . . . 

(bpa*caa+bpe*deea*ceaa) , bpe*cee, (ap-bpa*daa*cap. . . 

-bpe *deea * deaa * cap-bpe* deea *ceap-bpe*dee* cep . . . 

♦bpa * daya* cyap+bpe* deea* deaya *cyap+bpe* deye* cyep) ) ; 
btc- (baa; beea *deaa; (bpa *daa+bpe*deea* deaa) ) ; 
etc- ( ceaa, 0 ‘ones (pea, nke) , (-deaa* cap+deaya*cyap ) \ ; 
dtc- (deaa ) ; 


// Decompose the Input and output weighting functions into 
// their state space representations 
(awl, bwl, cwi,dwl) -split (swi,nwl) ; 

(awo, bwo, cwo, dwo) -split (swo, nwo) ; 


// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 


Compute the costs; fpO (performance), fpl (tracking) and their gradient 
fpO-Integra 1 (tr ((Wo(K- tilde (K) )W1) A T) (Wo (K-ti lde (K) ) Wi) ) ) 
f pi -I ntegra 1 (tr ((nrm(T- 1 1 lde (T) ) ) *T) (nrm (T-t 1 lde (T) ) ) ) 

DfpO-partit loning cost gradient 
— 2*Re( (LeftFactor) 

I " t tlnv(sl-A)*B; 1 1 *W1 * ( (Wo (K- 1 1 lde (K) ) Wl ) *T) *Wo‘ (C* lnv (s I-AI , I 
(Right Factor) )' 

Dfpl-tracklng cost gradient 
— 2*Re( (LeftFactor) 

Integral 



// ( ( inv (sI-A) *B; 1 1 * (nrm* (T- tilde(T))*T*nrm*(C*inv(sI-A) , I) ) 

// (RightFactor) )' 

// Initialize the costs and gradients to zero 

fpO-O; 

fpl-0; 

Df 0-0* ones (nka + nke + np+ma4me+ la + le f nka + nke + np + ka + ke) ; 

Df 1-0* ones (nka + nke + np+ma f nka + nke + np+pea) ; 

// Integrate over the frequency interval [FRO (1) , FRQ (2) ) using Simpsons rule, 
// coefficients are in the variable WEIGHTS.. 

// Build up the partitioning costs fpO, fpl, and the integral 
// part of the gradients DfO and Dfl 
// DfO-Integral 

// ( Ilnv(sl-A) *B; I ) *Wi * ( (Wo (K- ti lde (K) ) W l ) A T) *Wo* (C*inv (sI-A) , I) ) 

// Dfpl-Integral ( ( i nv (sI-A) *B; I)Mnrm*(T- tl lde (T) ) A T*nrm* (C* inv (sI-A) , I) 
for i-1 : f rq (3, 1 ) , om-omega (1) ; . . . 
kO-inv (om* jay*eye (nt ) -atO) ; ... 
wi-cwi*inv (om* Jay* eye (nwi) -awi) *bwi4dwi; . , . 
wo-cwo* inv (om* jay* eye (nwo) -a wo) *bwo4dwo; . . . 

kdif f -wo* (ck * i nv (om* Jay* eye (nk) -ak) *bk4dk-ct 0* k0*bt0-dt0) *wi ; . . . 
fpO-fpO 4welght (i) *om*sum(dlag (kdiff' *kdif f) /Nperf (1) ) ; . . . 

Df0-Df04 weight (1) *om/Nperf (1) * (k0*bt0;eye (mtot) | *wl*kdlf f' *wo*. . . 

[ct0*k0, eye (kt) ) ; . . . 
nrm-dlag ( (lambda . /Nt rack (i , : ) ) * * .5) ; . . . 
kl-inv (om* jay*eye (nt ) -ate) ; ... 

ktrc-nrm* (eg* Inv (om* jay* eye (ng) -ag) *bg4dg-ctc* kl *btc-dtc) ; . . . 
fpl- fpl 4 weight (1) *om* sum (dlag (kt rc' *kt rc) ) ; . . . 

Of 1 -Dfl 4 weight (i) *om* (kl *btc;eye (ma) ) *ktrc' *nrm* (ctc*kl, eye (pea) ) ; . . . 
end; 

A 

-'//It the airframe controller is to be optimized, calculate the cost gradients 
//with respect to the parameters belonging to the airframe controller by 
//multiplying DFO and Dfl by the corresponding Left and Right Factors 
//This is done for both the partitioning cost DfaO and the tracking 
//cost DfeO (with respect to the airframe parameters) 

//DfaO- ( (LeftFactor) DfO (RightFactor) )' 

//Dfal- ( (LeftFactor) Dfl (RightFactor) )' 
if aoreol, . . . 

DfaO- ( (eye (nka) , z (nka, nke4np4mtot) ; z (ma, nt) ,eye (ma) , z (ma,me4lt) ? . . . 
z (la, nt4mt ) , eye (la) , z(la,lc) |*DfO* (eye (nka) , z (nka, ka*pea) ; z (nke, nka4ka ) . . . 
beea ; z (np, nka) , bpa , bpe*deea; z (ka, nka) , eye (ka ) , z (ka, pea) ; . . . 
z(ke,nka4ka) , deea ))';.. 

Dfal- ( (eye (nka) , z (nka, nke4np4ma) ; z (ma, nka4nke) , -cap, eye (ma) ; z (la, nka 4 nke) , . . . 
cyap, z(la,ma) ) * D f 1 • (eye (nka) , z (nka, ka 4 pea) ;z (nke, nka4ka) ,beea; z (np, nka) , . . . 
bpa, bpe*deea; z (pea, nka4ka) , eye (pea) ))';.. 
else Df a0-0*ones (ska) ; Dfal-DfaO; end, 

//Similarly, if the engine controller is to be optimized, 

//calculate the cost gradients with respect to the engine 
//parameters by multiplying DfO and Dfl by the requisite Left and 
//Right Factors. Gradient of partitioning cost with respect to the 
//engine parameters is DfeO and gradient of the tracking cost is Dfel. 

//DfeO- ( (LeftFactor) DfO (RightFactor) ]' 

//Dfel- ( (LeftFactor) Dfl (RightFactor) )' 

1 f aore<>2, . . . 

DfeO- ((z (nke, nka) ,oye (nke) , z (nke, np4mtot) ?ceaa, z (pea, nke) ,-ceap, deaa, . . . 
z (pea, me) , deaya, z (pea, le) ; z (me, nt4ma) , eye (me) , z (me, It) ;z (le, nt4mt4la) , . . 
eye (le) 1 . . . 

*DfO*(z(nka,nke4ke);eye(nke),z(nke,ke);z(np,nke), bpe; z(ka,nke4ke) ;z (ke,nke) , . 
eye (ke) ))';... 

Dfel - (( z (nke, nka) , eye (nke) , z (nke, np4ma) ;ceaa, z (pea, nke) , ... 

- (deaa *cap4ceap-doaya * cyap) , deaa ; z (me, nka 4 nke) , -cep, z (me, ma) ; z ( le, nka 4 nke) , . 
cyop, z (le,ma) | *Dfl * (z (nka, nkefke) ;eye (nke) , z (nke, ke) ;z (np, nke) ,bpe; . . 
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z (pea, nke4ke) ))';... 
else Dfe0-0*ones (ske) ;Dfel-Df eO;end, 

//Assemble the final cost and gradients using the formulas 
//described earlier and the airframe and engine costs and gradients 
//Just generated. 

coef- (log (f rq (2, 1 ) / f rq (1 , 1) ) ) / (3* f rq (3, 1) ) ; 
fpO-real (coef *fp0) ; fpl-real (coef * f pi ) ; 
f p-fp04 fpl ; 

DfO-- 2* coef *rea 1 (long col (DfaO, DfeO) ) ; 

Df 1 — 2* coef * real (long col (Df al , Dfel ) ) ; 

Df=Df04Dfl; 

ret f 

// COST evaluates the total cost of approximating a centralized controller by 
// hierarchically partitioned subcontrollers. It Is the sum of two costs, 

// fPerf which measures the performance cost of approximation and fTrack which 
// measures the error in meeting the command tracking requirement. It's input 
// is the parameter vector, p, and the tracking weighting vector, lambda. 

// The file ' const.dat' contains the constants of the process 

// the dimensions ka, ma, ke, me, pea, np, nka, nke, nk, nwi (o) , ng, la, le (k-outputs, 

// m-inputs, 1- neg. feedbacks, n-order, p-intermediate variables) 

// the range vector of frequencies frq 

// the actual frequencies omega 

// the plant in system form sp 

// the global controller in system form sk 

// the input and output weighting matrices in system form swl, swo 

// the global tracking command matrix in system form sg 

// the weights used in the numerical Integration weight 

// the normalization matrix for the performance cost Nperf 

// the normalization matrix for the tracking cost — Ntrack 

// 

// The state space representation for the weighted difference between global 
// controller and assembled partitioned controller is determined. The 
// H_2 norm of the resulting system is computed. The norm is calculated by 
// applying a Simpson's rule suitably modified to account for the exponential 
// distribution of the omegas. 

// 

// The normalized H_2 norm of the difference between the tracking command 
// transfer matrix and the nominal one is also calculated. This is added to 
// the previous H_2 norm. The gradient of this sum is computed. 

// 

// The outputs are fpO the H_2 norm of the weighted difference between 

// centralized and partitioned controllers 

// fpl the normalized H_2 norm of the deviation from 

// nominal command tracking 

// fp - fpO 4 fpl 

// and DfO, Dfl and Df the gradients of fpl, fpl and fp resp. 


// Created by Phil Schmidt and students Nader Kamrani and Brian Holawecky at 
// The University of Akron with the support of NASA Lewis Research Center 
// under grant NAG-3-1146. 
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// (all, fal, fpO, fpl, dfdpal , dfpO, dfpl, flag) -INACCURATE (x, £z,dfdpz, . . . 
d, fmln, del ta f , lambda , a lphamax) 

// 

// The function INNACURATE In file LINESRCH . MTX 
// performs Fletchers Inaccurate line search as part 

// of the unconstrained optimization process. Given the search direction, 

// a point is found that satisfies the Arml jo-Goldsteln conditions using 
// a two stage bracketing/sectioning procedure. 

// Input: 

// x - (real,vec) the current point corresponding to al-0 

// fz - (real, seal) objective function value for al-0. 

H dfdpz - 1 rea 1 , voc) directional derivative of objective function at al-0 

// d - (real,vec) the direction vector which the BFGS procedure 

// has chosen to perform the univariate search. 

// fmln - (real, seal) user-defined lower bound on the problem; 

// if a point is found such that the objective value is less than 

// fmln, then the procedure will terminate. 

// deltaf - (real, seal) estimated change in cost 

// alphamax - (real, seal) maximum alpha set be stability constraint. 

// Output: 

// all - ( rea 1 , sea 1 1 the step-size determined to be the solution 
// to the univariate problem. 

// fa 1 , fpO, fpl - (real, seal) the cost function values at the point 
// ( x + a 1 * d ) ; this Is Included to help keep the number of computations 

// to a minimum. 

// dfdpal, dfpO, dfpl - (real,vec); gradients at the min point. 

II flag - ( lnt , seal ) denotes the condition with which INACCURATE has 
// completed the line search. 

// 0 - solution found; 1 - solution with function value less than fmln 

// 2 - reached end of search Interval, reduction without finding min 

// 


// The Armijo/Goldstein parameters, described in more detail in 
// Fletcher, PRACTICAL METHODS OF OPTIMIZATION, Wiley, 1987. 

Rho-0.01; 

Sigma -0.7; 

Taul-10; 

Tau2-0. 1; 

Tau3-0. 5; 

f pz-dfdpz' *d; // slope when alpha-0, ie directional slope 

a lpr-0; // alpha (1-1) the previous alpha 

faipr-fz; // function value at previous point 

fpalpr-fpz; // slope at previous point 

bounds- (l;- 10 *deltaf/ fpz; alphamax) ; // used to find inital alpha 
al-mln (bounds) ; // initial alpha 

albound»min((fmin-£z) / (Rho* fpz) , alphamax) ; // search restricted to (0,albound) 


// The following algorithm Is not documented in great detail as it follows 
// almost directly from Fletcher, PRACTICAL METHODS OF OPTIMIZATION, Wiley, 
// 1987 pp.34,35. The only part of the code not in the book is 
// the minimization of the cubic interpolation of fal, fpal, falp 
// and fpalp used to generate the next alpha in both the bracketing 
// and the sectioning phases 


// The bracketing phase 
whl le 1>0 do ... 

(dfdpal, dfpO, dfpl, fal, fpO, fpl ) -COST (x+al *d, lambda) ; ... 

fpal -d fdpa 1' *d; . . . 

al, fal, fpal, alpr, falpr, fpalpr,... 

if fal <-fmi n then flag-1; ali-al; retf; end;... 

if fal>f z+al *Rho* fpz, a -a lpr;b-al ; fa- falpr; fb-fal; fpa- fpalpr; fpb-fpal; . . . 
if X fai>faipr,’ a-alpr;b-a 1; fa-falpr ; fb-fa 1; fpa-fpalpr; fpb-fpa 1 ; exit ;end; . . . 
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If abs (fpa 1) o-Slgma* fpz, f lag-0;al l-al; retf ;end; . . . 

if fpa l>-0 , b-a 1 pr ; a -a 1 ; f b- f a 1 pr ; f a - f a 1 ; f pb- f pa 1 pr ; f pa - f pa 1 ; ex i t ; e nd ? ... 
if abs (a 1-albound) <0 .001 *albound, f lag-2; ali-al; retf ;end; . .. 
if albound<-2*al-alpr then m- (al+albound) /2; n-albound; . . . 
else m»2*al-alpr; n*min ( (albound, al+Taul* (al-alpr) ) ) ;end; . . . 

da l-al -alpr; dfal- (fal -falpr) /dal; d fpa 1 - ( fpal-f pa lpr) /dal;. . . 
c4- (dfpa 1-2* (dfal- fpa lpr) /dal) /dal; c3-0 . 5*d fpal-1 . 5*dal *c4; . . . 
if c4--0, if c3--0, alnew— (abs (fpalpr) /fpalpr) *albound; .. . 
else alnew- a lpr- fpa lpr/ (2*c3) ;end; . . . 

eisei f c3*c3-3*c4*fpalpr<-0, alnew-- (abs ( fpalpr) / fpalpr) *albound; . . . 
else alnew»al+(sqrt (c3*c3-3*c4*fpalpr)-c3)/ (3*c4) ; end;... 
alpr-al; falpr-fal; fpalpr-fpa 1 ; . . . 

If a lnewom, a lnew-m;elsei f alnew>-n, a 1 new-n; end, . . . 
al-alnew; . . . 
end; 

// The steps above select a new alpha in the interval 

// ( a 1 ♦ (a 1 -a lpr) ,al + Taul* (a 1 -a lpr) ) . alpha is chosen to minimize 

// a cubic polynomial which interpolates the values of f(al), f'(al)# 

// f (a lpr) and f' (alpr) . Note that the new point is chosen by moving far out 
// to the right in an effort to bracket an interval of acceptable points. 

// Similarly, in the sectioning phase, construct a cubic lnterpolant 
// using the values of f (A) , f (B) , f' (A) and f' (B) and find a minimum point 
// which (hopefully) lies in the interval 

// ( (1 -Tau2) # A+ Tau2*B,Tau3*AMl-Tau3) *B1, 0<Tau2<Tau3<l . The new 

// point is chosen in this interval (which is a subinterval of its 
// predecessor). Eventually an acceptable point will be found. 


// The sectioning phase 
while 1>0 do ... 
a, fa, fpa, b, fb, fpb, . . . 

da 1 -b-a; d fa 1- (fb-fa ) /da 1 ; dfpal- ( fpb-fpa) /dal ;.. . 
c4- (dfpa 1-2* (dial- fpa) /dal) /dal; c 3 - 0 . 5* df pal-1 . 5* da l* c4 ; ... 
if c4-=0, if c3--0 , alnew«b; else alnew-a-fpa/ (2*c3) ;end; . . . 
eisei f c3*c3-3*c4*fpa<-0, alnew-b;... 

else a lnew-a + (sqrt (c3*c3-3*c4*fpa)-c3) / (3*c4) ; end; . . • 
m-a+tau2*dal; n-b-t au3*da 1 ; . . . 

if m<n, if alnew<«m, alnew-m; elself alnew>-n, alnew-n; end,... 
else if alnewon, alnew-n; elseif alnew>-m, alnew-m; end, end;... 
(dfdpal, dfpO, dfpl, fal, fpO, fpl ) -COST (x + alnew*d, lambda); ... 
fpal-dfdpal' *d; . . . 

if fa 1 <-fmin , flag-1; ali-alnew; retf; end;... 
if fal>f z+Rho*alnew* fpz, b-alnew; fb-fa 1; fpb-fpal; .. . 
elself f a 1> fa, b-alnew; fb-fal; fpb-fpal; . . . 

else if abs (fpal) o-Sigma*fpz, flag-0; al 1-alnew; ret f; end;... 
aold-a; faold-fa; fpaold-fpa; a«alnew; fa-fal; fpa-fpal;... 
if 0<- (b-aold) *fpal, b-aold; fb-faold; fpb-fpaold; end;... 
end; . • . 
end; 


ret f 

// The function INACCURATE in the file ' 1 i nesrch .mtx' performs Fletcher's 
// Inaccurate line search as a part of the unconstrained optimization 
// procedure, PARTITION. It makes calls to the user-defined function (UDF) 
// COST in ' costO.mtx' which produces both the function value and its 
// gradient. 

// INACCURATE will solve the univariate minimization problem required 
// at each line search in the overall optimization process, that is, 

// 

// min f (x+a 1 *d) 

// This is accomplished by use of a two-stage (bracket ing/ sect lonl ng) 
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procedure which employs cubic approximation that will seek 

only to find a point which satisfies the Armi jo-Goldsteln conditions. 

This approach is useful because the BFGS method of 
determining a search direction has proven to be robust enough to 
be effective even with an Inaccurate line search. 

This line search makes use of many parameters that define various 
aspects of the search computations. These parameters have been set 
to values that have proven to be adequate for this procedure. Of 
course, these may be changed to suit the user's needs. For a more 
thorough description of how each parameter is used in this procedure, 
see Fletcher, PRACTICAL METHODS OF OPTIMIZATION, Wiley, 1987. 

Parameters: The Armi jo/Goldstel n parameters should satisfy: 

0 < Rho < 1/2 smaller is easier to satisfy; we use 0.01 

Rho < Sigma < 1 Larger is easier to satisfy; we use .7 

The next is used in the bracketing stage to find (if 
necessary) larger and larger intervals which might contain an 
acceptable, i.e. Armi jo/Goldsteln, interval: 

Taul > 1 we use 10, the larger, the more effort might be 

used in the sectioning phase. 

The last two parameters are used to reduce the acceptable 
interval until an acceptable point is found: 

0< Tau2 < Sigma is advised wo use 0.1 

Tau2 < Tau3 <- 1/2 we use .5 to get the greatest 

reduction of size. 

Notes on the use of the output variable 'flag' in this function: 

Flag is used by INACCURATE . DAT to characterize the condition of 
the line search when it terminates. The values can signify: 
flag-0 - the solution to the line search has been found; 

nothing special to note about the solution. 
flag-1 - a point has been found that gives an objective value 
less than fmln. 

flag-2 - the bracketing step has converged to the right hand 
endpoint of the Interval without satisfying the A-G 
conditions. It should be noted however that the new point 
will have a reduced function value (although, perhaps 
insufficient to trigger the A-G conditions). 


// Created: 01/30/90 

// Programmer: Steven Ims 

// Rewritten by Phil Schmidt and students Nader Kamranl and Brian Holawecky at 
// The University of Akron with the support of NASA Lewis Research Center 
// under grant NAG-3-1146. 
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// ( sm J -modi (s, ns) 

// The function MODL In file MODL.MTX puts a system matrix Into a modified 
// modal form reducing the number of parameters In the optimization. 

// The Inputs are a system matrix and Its order. The system Is put Into the 

// form where A has 2x2 companion matrix blocks whose first rows are (0 1) and 

// whose second rows are [a b] . The transformation matrix Is normalized by 
// requiring that all nonzero entries In the first column of the B matrix 
// remain fixed (In case of zero entries In Bl, the entry Is fixed at le-9) . 

// The C, and D are full. This should work except in the special case of 

// repeated real roots. 

// Input: 

// s - the system matrix to be put In modal form 
// ns - the order of the Input matrix 
// Output : 

// sm - the system matrix, now In modified modal form, still has order ns 
// 


// Convert the input system matrix S into the internal Matrixx modal form 
// where the A submatrix Is of the form that the real eigenvalues 
// are on the diagonal and the complex conjugate eigenvalues (a + bl) 

// and (a - bl) are stored in 2x2 matrices of the form (a b) on the first 
// row and [-b aj on the second row. 

( sm, t ) -modal (s, ns) ; 


// Split the modal form sm into Its components 
(am, bm, cm, dm) -s pi It (sm, ns) ; 


// Create a permutation matrix P which collects all the real distinct 
// eigenvalues at the top of am and sends all the 2x2 matrices representing 
J\ // the complex eigenvalues to the bottom of am. 

! ^ 1-1; j-0; p-eye(ns); 

( n, o) -size (cm) ; 
while 1 < ns, . . . 
if am(i, 1+1) -0, 1-1+1;... 
else, . . . 

m-p ((1:1 + 1),:);... 

If l>j + l, for k-1 :-l: j+2, p (k + 1 , : ) -p (k-1, : ) ? end; end;... 

P((j+1: j+2),:)-m;... 
j-j+2; 1-1+1;... 
end; . . . 
end; 


// Rearrange the real eigenvalues so that the ones that are close 

// in value are separated. 

nr-ns-j; 

If nr > 2, q-mod(nr,2); nrhal f- (nr-q) /2; tpnt-J + 2; mpnt-J+nrhalf+1 
for 1-1: (nrhalf-1) , temp-p (mpnt, : ) ; . . . 

p( ( (tpnt + 1) :mpnt ) , :) -p((tpnt : (mpnt-1) ),:);...* 
p (tpnt , : ) -temp; tpnt-tpnt+2; mpnt-mpnt +1; . . . 


end; . . . 
end; 


// Since S-(A,B; C,D) is the form of a system matrix, multiplying 

// (P, i ) * S * (inv(P); I| will have the following effect on A,B,C,D 

a-p*am* inv (p) ; 

b-p*bm; 

c-cm* lnv (p) ; 

ft If two real eigenvalues are still very close in value, shift one slightly 
// more negative 
for i-l:nr-l, ... 

If a (1, 1) -a (1 + 1, 1 + 1) , a (1, 1) -a (1, 1) -10e-9; end; 
end; 
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// Construct a transformation matrix T such that 

// T * a * inv (T) has the effect of taking 2x2 blocks of 

// real eigenvalues (or of complex conjugate pairs) and 

// returning the desired structure (0 1; a b) which have the same 

// eigenvalues as the 2x2 blocks. 

t-0*eye (ns) ; 

q-mod (ns, 2) ; 

for 1-1 : (ns-q) /2, . . . 

t (2*1-1, 2*1-1) -1; t (2*1-1, 2*1) -1; ... 
t ( 2 * 1 , 2 * 1 - 1 ) -a ( 2 * 1 - 1 , 2*1-1) -a (2*1-1, 2*1);... 
t (2*1, 2*1) -a (2*1, 2*1) -a (2*1, 2*1-1) ; ... 
end; 


// The variable q keeps track of an odd size matrix, i.e. an odd entry which 
// does not fit In the 2x2 blocks 
1 f q-1, t (ns, ns) -1; end; 

H (T 1) * system matrix * (lnv(T); I) will generate the following 
t lnv-inv (t ) ; 

aa-t*a*tinv; bb-t*b; cc-c*tlnv; 


// Set the first entry In the modal blocks to zero; (0,l;a.b) 
for 1-1 : (ns-q) /2, aa (2* 1-1, 2* 1-1) -0; end; 

// Construct the system matrix sm which is now In the modified modal form 
sm- (aa, bb; cc, dm) ; 


ret f 

// Created by Phil Schmidt and students Nader Kamranl and Brian Holawecky at 
// The University of Akron with the support of NASA Lewis Research Center 
// under grant NAG-3-1146. 



// [ska, ske) -mat (p) 

// 

// The function MAT In file PARMAT.MTX generates 
// the partitioned system matrices ska and ske from 

// the parameter vector p. Note that p Is: Aa, Ae, Ba,Ca, Deaa,Deaya,Be,Ce,Deea. 
// The other submatrices are all constants and are loaded from const.dat 
// Input: 

// p - the vector of parameters 
// const.dat - datafile of constants 
// Output: 

// ska - the airframe controller 
// ske - the engine controller 
// 

// load the default data from const.dat 

load 'const.dat' nka ma lya nke bal ka pea me ke lye bel ... 
daa daya dee deye skaO skeO aore fixd 

// Set the counter, first to 0 then to nka etc. as new blocks are built from p 
dum - 0; 

// If airframe is not fixed, then create a blank aa (A submatrix for ska) 

// and copy the requisite entries from the parameter vector p. Recall the 
// modal from of the system matrices (refer to modl.mtx) reduces the 
// A submatrices to 2x2 companion blocks of the form (0,1; a,b), so only 
// the a and b need to be read 
if aoreol, aa-0*eye (nka) ; q-mod (nka, 2) ; ... 
for i-1 : (nka-q) /2, . . . 

aa (2*1-1, 2*1) -1; aa (2M , 2*1-1 ) -p (2* 1-1, 1) ; ... 
aa(2*i,2*i)-p(2*i,l); end;... 
jy if q-1, aa (nka, nka) -p(nka, 1) ; end; ... 
j\ matot-ma+lya; katot-ka+pea; dum-nka;... 
end, 

// If engine is not fixed, then create a blank ae and copy requisite 
ft entries from the p vector, with the same provisos as for aa above 
if aore<>2, ae-0*eye (nke) ; q-mod (nke, 2) ; ... 
for i-1 : (nke-q) /2, . . . 

ae (2*1-1, 2*i) -1; ae (2M , 2*i-l) -p (dunU 2* 1-1, 1) ; ... 
ae(2M,2*i)-p (dum+2*i, 1) / end; . . . 
if q-1# ae (nke, nke) -p (dum+nke, 1) ; end; ... 
metot-me+pea+lye; dum-dum+nke; . . . 
end, 

// If airframe is not fixed, generate the ca,deac,da entries from p. 

// Note that the first column of the ba entry is fixed as required by 
// the modal form, so load it from bal (saved in const.dat) 

// The rest of ba,ca,deac are loaded from p. If Ds are not fixed, then 
// their values are read from p, otherwise from constants 
// (and daya depending on whether feedback lya exits) . 
if aoreol# ba (: , 1) -bal; . . 

for i-2:matot, ba (:, i) «p ( (dum+ 1 : dum+nka ) ,1) ; dum-dum+nka; end;., 
for i-l:nka, ca ( : , 1 ) -p { (dum* 1 :dumf katot ) , 1 ) / dum-dumf katot ; end/., 
if (fixd-1) * (fixd-2) <> 0, ... 

for i-l:matot, daaya (: , i) -p( (dumH :dum*ka) , 1) ; dum-dum+ka; end;., 
else if lya-0, daaya-daa; else daaya-[daa, daya); end;... 
end; . . . 

for i-l:matot, deac(:, 1) -p ( [dum+l :dum+pea) , 1) ; dum-dumfpea; end;., 
da- (daaya; deac) ; . . . 
end, 

// If the engine controller is not fixed, load the first column of be 
// from the bel entry stored in const.dat. Then load rest of be,ce,deea 
// from p. Construct de from deea, dee (and deye depending on whether 


parmat.mtx 

// there is any feedback lye) . 
if aore<>2, be ( : # 1 ) -bel ; . . 

for i-2:metot, be (:, i) -p ( (dum*l:dum+nke) , 1) ; dum-dum+nke; end;., 
for l-l:nke, ce (:# i) «p ( [dum+1 :dum4ke) , 1) ; dum-dum+ke; end;., 
if (f ixd-2) * (f ixd-3) -0, deea«0*ones (ke, pea) ; . . . 

else for i-l:pea, deea ( : , i ) -p( (dum+l :dum+ke) , 1) ; dum-dum+ke; end;... 
end; ... 

if (f ixd-1) * (fixd-2) <> 0,meplye-me4lye; . . . 

for i-l:meplye, deeye (:, i) -p ( [dum+1 :dum+ke) , 1) ; dum-dum+ke; end;., 
else if lye-0, deeye-dee; else deeye-(dee, deye); end;... 
end; . . . 

de-[deea, deeye);... 
end, 

// If the airframe controller is not fixed, then construct ska from aa,ba 

// ca and da. If there is no feedback defined lya-0, then extend 

// ska by one column of zeros to accomodate the dummy feedback that will 

// be used. 

// If the airframe controller is fixed, then use ska-skaO 
if aoreol, ska- [aa, ba; ca, da ) ; . . . 

if lya-0, (rska, cska J -si ze (ska) ; ska- [ ska, 0*ones (rska, 1 ) ) /end, . . . 
else ska-skaO;... 
end, 

// If the engine controller is not fixed, then construct ske. 

// If there is no feedback defined lye-0, then add the zero column 
// to accomodate the dummy feedback. 

// If the engine controller is fixed, then let ske-skeO the fixed 

// engine controller. 

if aore<>2, ske- (ae, be; ce, de) ; . . . 

if lye-0, ( rske, cskej -si ze ( ske) ; ske- (ske, 0*ones (rske, 1 )] ; end,... 
else ske-skeO; end, 

ret f 

// This program takes the parameter vector P, which is a column vector, and 
// uses the dimensions of the subcontrollers stored in 'const.dat' to 
// reconstruct the partitioned system matrices ska and ske. It allows for the 
// possibility that one of the subcontrollers is fixed and loads its initial 
// value. Note that the constant submatrices Bal, Bel, Daa, Daya, Deac, Dee and 
// Deye are properly loaded. 

// Created by Phil Schmidt and students Nader Kamrani and Brian Holawecky at 
// The University of Akron with the support of NASA Lewis Research Center 
// under grant NAG-3-1146. 
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// p_i-partitlon (stop) 

// The function PARTITION In file PARTITIO.MTX Is the 
// outer level loop in the optimization process. 

// It Implements the Broyden-Fletcher-Goldfarb-Shanno method of finding 
// a search direction for ml nlml zat Ion, then calls function INACCURATE 
// to generate a new point by minimizing in that direction. This 
// process is then repeated until the convergence conditions (as 
// checked by function CONVERGE) are met or some other stopping criteria 
// are met (number of iterations exceeds maximum or function value is 
// less than fmin) . 

// Input: 

// stop - the vector of stopping conditions 
// const.dat - the datafile of constants 

// inter.dat - the datafile containing intermediate results 

// Output: 

// pi - the ith point along the optimization process 

// inter.dat - the datafile where intermediate results are stored 

II 


load 'const.dat' omega nka nke atot stabil 

load 'inter.dat' p i gradl gradO gradl fh JhO j hi hi lambda fx; 

( row, col ) -size ( p i ) ; (cnt , tmp) -si ze (f h) ; 

// Prestore the messages to save the length of loops, in order that the 

// Matrixx line limit be avoided 

Is-' LINES EARCH TERMINATED BY HITTING BOUND'; 

not-' NOTICE! STABILITY BOUND.'; 

a_ph-'A phat has been found that gives an objective value < fmin.'; 
ex-'exceeded maximum number of Iterations'; 
plpar-' STRIP YLAB/f (P) I fO (P) I f 1 (P ) /' i 

// Read stopping conditions from the vector stop 

epsl -stop(l) ‘ones (p_i) ; 

del -stop (2); 

eta - stop(3); 

maxcnt -stop(4); 

fmin -stop(5); 

dfx-dlag (ones (row, 1) - ( fx; 0*ones ( (row- a tot) ,1) ) ) ; 
del f -max (0. 9*abs (fh (cnt) ),10**(-6)); 
cvgfact-1 ; 
cvf lg-0; 

// While the convergence/stopping conditions have not been met 
while cvflg<0.5 do... 
cnt, fh (cnt) , ... 
if cnt>l, . . . 

if cnt>20, ITER- ( cnt-20 : cnt ) ' ; ... 
else ITER- ( 1 : cnt J ' ; . . . 
end . . 

Y- ( fh (ITER) , jhO (ITER) , jhl (ITER) ) ; PLOT (ITER, Y, plpa r) ; . . . 

end, ... . .. 

di— dfx*hi *dfx*gradi; . . . // search direction is -Approx Inv Hessian gradient 

alphamax-lel 6; . . .// check on stability constraints 
for 1-1 : atot , . . . 
if fx (I ) *0, ... 
if dl ( I) >0, . . . 

if (stabl l-p_i (I) ) /di (I) <alphamax, . . . 

alphamax-(stabil-p_l (I)) /di (I) ; . •• 
end, . . . 
end, . . . 

else di (I)— 0; end,... 

!all*fp_ll,J0, jl, gradil, gradO, gradl, flag) -INACCURATE <p_l,fh (cnt), gradl, ... 


.mtx mmm 

dl, fmin, del f, lambda, alphamax) // calls the llnesearch 
if flag-2, display (Is) ; end;... 

if ali-alphamax, cvgfact-0; di splay (not ) ; end,... 
fh (cnt +1 ) -fp_i 1 ; . . . 
p_il-p_i+ali*di; ... 
jhO (cnt + 1) -j0; ... 
jhl (cnt+1) - jl; . . . 

cvflg-CONVERGE(p_ll, fh (cnt + 1 ) , p_i , f h (cnt ), gradil , cvgfact *epsl , . . . 

cvgfact *del, eta, fx) ;...// checks convergence 
cvgfact-1; . . . 
psi-ali*di;... 

Qi-dfx* (gradil-gradi ) ; 
psiQl-psi' *Qi; . . . 

pQmat- (eye (row) -psi *01 ' /pslQi) ; . . . 

hi -pQmat *hl * (pQmat ' ) + ( (psl # psl ' ) /pslQl) ; . . .// updates Inv Hessian Approx 
for I-l:atot, ... 

if fx (I ) -0, if p_il (I)>-1 .001*stabll, fx (I) -1; . . . 
p 11 (I)-stabll;dfx (I, I) =0; end. . . 

elseif gradil (I)>0, fx (I) -0; for J-l : row, hi (I, J) -0;hi (J, I) -0;end, . . . 

hi (I, I) -l;dfx (I, I) -1; . . . 
end, end ...II updates stability constraints 
p_i-p_ 11 ; . . . 
gradi-gradil ; . . . 

del f -max (fh (cnt) -fh (cnt + 1) , 10* * (-6) ) ; . . . 
cnt-cnt+1 ; . . . 
if mod (cnt , 20) -0, ... 

save 'inter.dat' p_l fh jhO jhl gradi gradO gradl lambda fx hi;end... 
if cnt -maxcnt, disp (ex) ;.. . 

save 'inter.dat' p_i fh jhO jhl gradi gradO gradl lambda fx hi;... 
ret f ; end, . . . 
if f lag-1, disp (a_ph) ? 
save 'inter.dat' p_l 
ret f ; end; . . . 

save 'inter.dat' cnt p_i fh jhO jhl gradl gradO gradl lambda fx hi; 
d i sp (' Convergence of f-value and parameters has occurred'); 
if cvflg-2, dlsp('all partials are also less than'); eta, end; 
etf 

'/ PARTITION is the top-level function in the optimization routine. 

7 It makes calls to the user-defined functions INACCURATE and CONVERGE 
7 in the files ' 1 incsrch.mt x' and ' converge .mtx' resp. 

7 It also uses the MATRIXX PLOT function to show intermediate results. 

t / 

7 PARTITIO.MTX is an implementation of the Broyden-Flet cher-Goldf arb-Shanno 
7 method of determining a search direction, with Fletcher's Inaccurate 
7 line search being used to locate the the next point in UDF INACCURATE. 

7 

7 PARTITION .MTX will solve hierarchical partitioning problems. 

7 

II Inputs: 

II stop * ( rea 1 , vec) column vector of stopping conditions 

II stop - (epsl del eta maxcnt fmin) 


column VKLIUL C i c * I V. — 

limits for elements of lp_ll - p_i I ; 
convergence of objective values, 

If (p_il) - f (p_i) I; 
convergence of gradients, 

| grad f (p_il) I < eta; 

maximum number of iterations of procedure 


fh jhO jhl gradi gradO gradl lambda fx hi; 


II 

Pa rameters 

being defined 

II 

START 


// 

II 

epsl 

(defaul t-le-9) 

II 

// 

del 

(default-le-9) 

// 

// 

eta 

(defaul t-le-9) 

// 

// 

maxcnt 

(default-100) 



// 

// 

// 

// 

// 

// 

// 

// 

II 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 

// 


allowed In searching for a solution, 
fmln (default-0) - lower bound on the objective function. 

Any point, x, found with an objective 
value less than fmln will be considered 
a solution. Since we minimize a norm, the 
default Is zero. 

Notes on the use of 'flag' in this function: 

Flag Is used by the INNACURATE functions to 

characterize the objective function and the current point, x. 

It can take on values as follows: 

(1) flag-0 - a point has been found which satisfies the Armijo- 

Goldstein conditions (see ' 1 inesrch.mtx' ) ; the precedure 
cont inues. 

(2) flag-1 - a point has been found that gives an objective value 

less than fmln. Procedure terminates. 

(3) flag-2 - a point has been returned by INACCURATE which doesn't 

satisfy the A-G condition, but which corresponds to the 
maximum allowable step In the search direction. This can 
happen because the steps are constrained to maintain 
subcontroller stability. Currently, the code continues 
as though the A-G conditions were satisfied. If problems 
such as looping are encountered, the user might try to 
reset HI to the identltiy at this point. 

The program terminates properly under one of the following conditions: 

(1) flag-1 - the objective value is less than fmln. 

(2) cnt-mxcnt - the maximum number of Iterations has occurred 

(3) cvflg-1 - the max deviation In parameters is less than epsl and 

the max deviation in COST is less than delta 

(4) cvflg-2 - same conditions as cvflag-1 and max deviation In 

partials is less than eta. 

On termination (of any type) the following can be found in the file 
' inter.dat' : 

p_i - last parameter vector (can use procedure MAT in ' parmatO.mtx' 
to generate the corresponding SKA and SKE) 
fh - complete total cost history (all Iterations) 

jhO, jhl - complete cost histories of fPerf and fTrack costs 
gradl - gradient of total cost at last parameter point. 
gradO, gradl - gradients of fPerf and fTrack at last point 

hi - last inverse Hessian updated according to values of p_i, etc. 
lambda - lambda vector 

fx - vector Indicating which entries In Aa and Ae are fixed at 
stability bounds 


// Created: 01/23/90 as BFGS.DAT 

// Programmer: Steven Ims 

// Revisions by Phil Schmidt and students Nader Kamranl and Brian Holawecky at 
// The University of Akron with the support of NASA Lewis Research Center 
// under grant NAG-3-1146. 
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// p-longcol (ska, ske) 

// 

// The function LONGCOL In file PARVEC.MTX 

// creates the long column vector of parameters p from the partitioned 
// system matrices ska and ske 
// Input: 

// ska - the airframe controller state-space matrix 
// ske - the engine controller state-space matrix 
// const.dat - datafile of constants 
// Output: 

// p - the long column vector of the parameters 

// 


load 'const.dat' nka nke me pea lye ma lya ka aore flxd 


(aa , ba, ca,da) -spl It (ska, nka) ? 
[ ae, be, ce,de) -split (ske, nke) ; 

metot»me+pea+lye; 
matot*ma+ lya? 

p-(0) ; 


// split ska and ske Into their 
// component (A,B;C,D) 

// define sires of the parameters 
// initialize p 


// If airframe controller is not fixed then copy entries that need to be 
// optimized from aa. Note that because aa has the modified modal from 
// consisting of 2x2 companion matrices (0,1 ? a,b), only the a,b need to 
// be included Into the parameter vector. The last block of aa may contain 
// only la,b) and needs to be Included If a Is odd (q-1) 
if aoreol,q“mod(nka, 2) ? . . . 

for 1-1: (nka-q) /2, p-(p;aa (2*1, 2*1-1) ? aa(2M,2 1)); end;... 

/\ if q-1, p- (p?aa (nka, nka) ) ? end?.. 

O end, 

// If engine controller Is not fixed, put the required elements of 
// ae into p (same modified modal form as aa) 
if aore<>2, q-mod (nke, 2) ? . . . 

for 1-1 :(nke-q)/2,p-lP;ae(2*l, 2*1-1) ; ae 12*1, 2*1)1 ;end;... 

If q-1, p-(p;ae (nke, nke) ) ? end?... 
end, 

// If airframe controller Is not fixed, then copy ba (note first column 
// of ba Is fixed and Is not copied), ca and da Into p 
i f aoreol, . . . 

for l-2:matot, p-(p;ba (: , 1) ) ? end?.. 

for 1-1 : nka, p-(p?ca (: , 1) 1 ? end?.. ...... ... .. 

if (flxd-1) * (f lxd-2) <> 0, for l-l:matot, p- (p?da ((1 :kal , 1) j , end, end,... 

for 1-1 :matot , p- (p; da ( (ka+1 :ka*pea ) , 1 ) ) ? end?... 
end , 

// If the engine controller Is not fixed copy be (first col fixed), ce and 
// de into p 
If aore<>2, . . . 

for l-2:metot, p-(p;be (: , 1) ) ? end?... 
for i-1 : nke, p-(p?ce (: , 1) ) ? end?... 

If ( f lxd-2) * ( f lxd-3) <> 0, for l-l:pea, p- (p?de {: , 1) 1 ? end, end,... 

If (f lxd-1) * (f lxd-2) <> 0, for 1 -pea+1 :metot , p-(p;de ( : , 1) ) ? end; end?... 

end, 

// Get rid of the 0 which was the 1st entry In p (used to Initialize p) 

(x y ) -size (p) ? 
p-p ( [ 2 : x ) , 1) ? 


ret f 


parvec.mtx 



// This program accepts the partitioned system matrices of subcontrollers 
// and generates the long column vector of the parameters. The order with 
// which this procedure builds p Is: Aa, Ae, Ba, Ca, Deaa, Deaya, Be, Ce, Deea . Note 
// that it skips matrices if the subcontroller Is fixed. It also Ignores the 
// submatrices which are constants Bal, Bel, Daa, Daya, Dee, and Deye. 


// Created by Phil Schmidt and students Nader Kamranl and Brian Holawecky at 
// The University of Akron with the support of NASA Lewis Research Center 
// under grant NAG-3-1146. 




// (ska_opt, ske_opt ) -restart (lambdanew, stop, newh) 

// 

// The function RESTART in file RESTART. MTX is used to 

// restart the optimization process after it has been intentionally 

// halted or has stopped due to the stopping criteria. 

// Some parameters may be altered if desired. 

// Input: 

// lambdanew - a new value for lambda, the weighting of the tracking 

// cost as part of the total cost function 

// stop - (optional) the vector of stopping conditions 

// newh - (optional) if any parameter whatsoever is passed, this sets 

// the inverse Hessian approxlmat ion to an identity matrix 

// const.dat - the datafile of constants 

// inter.dat - the intermediate results at the last checkpoint or the 

// final results if convergence/stopping criteria were met. 

// Output: 

// ska_opt - the optimized airframe controller 

// ske_opt - the optimized engine controller 

resize (' sstack' , 1000000) ; // load the relevant data 

load 'const.dat' nka nke pea lya lye; atot-nka*nke; 

load 'inter.dat' p_l fh jhO j hi gradi gradO gradl hi fx lambda; 

short e; // Specify the format for displaying matrices 

// Define UDF INACCURATE which implements Fletchers inaccurate line search 
define ' 1 inesrch.mtx' ; 

// Define UDF CONVERGE which checks if convergence conditions are met 
define ' converge. mtx' ; 

J\ 

O // Define UDF PARTITION which is the outermost loop and implements the 

// Broyden-Fletcher-Goldf arb-Shanno technique for finding a search direction 
// for optimization 
define ' parti tlo. mtx' ; 

// Define UDF Z which returns a zero matrix of the desired size 
define 'zero. mtx' 

// Define the UDF MAT which generates the SKA and SKE subcontrollers 
// from the parameter vector p 
define 'parmat.mtx' 

// Define UDF LONGCOL which generates the parameter vector from the 
// SKA and SKE groupings of the parameters 
define 'parvec.mtx' 

// Define UDF COST which evaluates the performance cost, tracking cost 
// and their sum, as well as the respective gradients for a specified 
// parameter vector 
define 'cost. mtx' 

// if stop does not exist, load it from par.dat 
if exist (' stop' ) -0, load 'par.dat' stop 

// Verify that the newlamda provided is a row vector or a scalar and is 
// consistant with pea 
( r lambda, c lambda ) -size ( lambdanew) ; 

if r lambda > 1, dlsp(' ERROR: lambda must be a row vector or a scalar!'), retf 
if clambda-1, lambdanew-lambdanew*ones (1, pea) ; ... 
elsclf clambda <> pea,... 

display (' ERROR: lambda must contain pea entries!'); retf; 
change 1 am-norm ( lambda- lambdanew) ; 
lambda -lambda new; 


restart. mtx 


// If the difference between old lambda and new lambda is large enough, 

// >le-10, then call COST to initialize costs and gradients, otherwise 

// just use old data 

if abs (changelam) >le-10, ... 

(gradi, gradO, gradl, fh, jhO, j hi ) -COST (pi, lambda ) ; 
else ... 

(cnt, tmp) -size (fh) ; f h (cnt) - jhO (cnt) ♦ jhl (cnt) ; gradi-gradO+gradl; . . . 
end, 

// If newh is defined, then set the Inverse Hessian to an identity matrix, 

// If any of the A parameters are at the bounds (noted in fx) 

// then set that element of hi to be zero (no further decrease in that 
// direction/parameter) 

if 1-exist (' newh' ) , ( row, col ) -size (p_l) ; hl-eye(row); ... 

for i«l:atot, if fx(l)-l, hi ( 1 , i ) —0 ; end, end,... 
end; 

// save the possibly new stopping vector and save intermediate data 
save 'par.dat' stop; 

save 'inter.dat' p_l fh jhO jhl gradi gradO gradl hi fx lambda; 

// call the PARTITION routine which eventually returns the final 
// optimized vector 
p_opt -PARTITION (Stop) ; 

// Generate the optimized ska and ske as the final output of the program 
(ska_opt, ske_opt ) -MAT (p_opt) ; 

load 'inter.dat' p_l fh jhO jhl gradl gradO gradl hi fx lambda; 

// If a dummy feedback for ya or ye was added, then remove it 
if lya-0, ( r, c) -size (ska opt) ; ska_opt-ska_opt ( : , 1 : c-1) ; end; 
if lye-0, (r, c) -size (ske opt) ; ske_opt-ske_opt ( : , 1 : c-1) ; end; 

save 'inter.dat' p_l fh jhO jhl gradl gradO gradl hi fx lambda ska_opt ske_opt 

// INTER.DAT now contains the final results including the optimal 
// subcontrollers. 

retf 

// Created by Phil Schmidt and students Nader Kamranl and Brian Holawecky at 
// The University of Akron with the support of NASA Lewis Research Center 
// under grant NAG-3-1146. 


// [skaopt, ske_opt J -start (lambda, stop) 

// The function START in file START. MTX Initiates the 
// controller partitioning optimization code. 

// It sets up initial data and precalculates some required information, 

// calls COST to initialize the costs and gradients and 
// calls PARTITION to do the optimization. 

// I nput : 

// lambda - a scalar or pea x 1 vector of weightings for determining the 

// relative contribution of fl(p), the tracking cost, to the total cost. 

// stop - (optional) the vector of stopping conditions 
// stop-(eps;delta;eta; iter; fmln) where 

// default: stop- (le-7; le-7; le-7; le2; le-2) ; 

// eps - tolerance for max change in parameters giving convergence 

ll delta - tolerance for max change in cost giving convergence 

// eta - tolerance for max change in norm of gradient giving convergence 

// iter - number of iterations to run; defaults to 100 

// fmln - minimum tolerance for the total cost; a cost below this 

// is assumed to be a minimum. 

// INIT.DAT - the data file containing the initial information 
// Output: 

// ska_opt - the optimized airframe subcontroller 

// ske_opt - the optimized engine subcontoller 

// 

resize (' sstack' , 1000000) ; // System stack size is Increased 

load * Init.dat" ; // Initial data is loaded 

// Define UDF MODI, which puts a system matrix into a modified modal form 
define "modi. mtx"; 

// Define UDF INACCURATE which implements Fletchers inaccurate line search 
define " llnesrch.mtx’ ; 

// Define UDF CONVERGE which checks if convergence conditions are met 
define " converge .mtx' ; 

// Define UDF PARTITION which is the outermost loop and implements the 
// Broyden-Flet cher-Gold farb-Shan no technique for finding a search direction 
// for optimization 
define ' part itio. mtx' ; 

// Define UDF RESTART which can restart the program using checkpointed 
// information 
define " restart .mtx' ; 

// Define UDF Z which returns a zero matrix of the desired size 
define "zero. mtx' 

// Define the UDF MAT which generates the SKA and SKE collections of the 
// parameters from the parameter vector 
define "parmat.mtx" 

// Define UDF LONGCOL which generates the parameter vector from the 
// SKA and SKE system controllers. 

Define "parvec.mtx" 

// Define UDF COST which evaluates the performance cost, tracking cost 
// and their sum, as well as the respective gradients for a specified 
// parameter vector 
define "cost. mtx" 

short e; // Specify the format for displaying matrices 

ska opt “0; ske_opt-0; // Initializing answers in case of an abrupt error/end 


if exist ("sp" )-0, i f exist (" spaug" ) -1, sp-spaug; 
if exist (" np" ) -0, if exist (" nspaug" ) -1, np-nspaug; 

// Checking the input data to see if it is consistant. 

if exist (" s ka")-0, display (" ERROR: s_ka is missing from init.dat"); retf ; 
if exist ("ns ka " ) -0, dl splay (" ERROR: nsjca is missing from init.dat"); retf; 
if exist (" s ke" ) -0, display (" ERROR: s_ke is missing from init.dat"); retf; 
if exist ("ns_ke')-0, display (" ERROR: ns_ke is missing from init.dat"); retf; 
if exist (" sc" )-0, display (" ERROR: sc is missing from init .dat "); ret f; 
if exist ("nsc" ) -0, display (" ERROR: nsc is missing from init .dat "); ret f; 
if exist (' sp" ) =0, display (" ERROR: sp is missing from init .dat "); ret f; 
if exist (" np" ) -0, display (" ERROR: np is missing from init .dat "); ret f; 
if exist (" frq" ) -0, dl splay (" ERROR : frq is missing from init .dat "); ret f; 
if exist ("pea ' ) -0, di splay (" ERROR : pea is missing from init .dat "); ret f; 

If exist i* la" ) -0, lya-0; else lya-la; 
if exist (" le" ) -0, lye-0; else lye-le; 

(r, c) -size (ska) ; ka_tmp - r-ns_ka-pea; ma_tmp - c-ns_ka-lya; 

[ r, c) -s 1 ze (s_ke) ; ke_tmp - r-ns_ke; metmp - c-ns_ke-pea-lye; 

(r, c) -size (sp) ; 

if ronp+ma__tmp+me_tmp*pea + lya + lye, ... 
display (" ERROR: number of rows in sp not consistant with other data ); ... 

retf; end; 

if conp+ka_tmp*ke_tmp, ... 

display (" ERROR: number of columns in sp not consistant with other data ),... 

retf; end; 

(r, c) -size (sc) ; 

if ronsc + ka_tmp+ke_tmp, ... 

display (" ERROR: number of rows in sc is not consistant with other data ); . 
retf; end; 

if consc+ma_tmp+me_tmp* lya + lye, ... 

dl splay (" ERROR: number of columns in sc not consistant with other data'); . 
retf; end; 

if frq (3,1) <- 0, ... _ . , . 

dl splay (" ERROR: number of observation points frq(3,l) must be positive ), . 

retf; end; 

if mod (frq (3, 1) , 2) -0, ... 

display (" ERROR: number of observation points frq(3,l) must be odd ); ... 

retf; end; 

if f rq (1, 1) >- frq (2,1), ... .. 

display (" ERROR: min. observation point frq(l,l) must be less than the ),.. 
di splay (" max . observation point frq(2,l)"); retf; end; 


// Get row and column of the centralized controller sc 
Irsc, esc) -size (sc) ; 

// The input and output weighting matrices swi and swo of orders nwi and nwo 
// resp. can be defined in init.dat. If not defined, they are assumed to 
// be identity matrices of the desired size (depending on the size of 
// the centralized controller) 

if exist (" swi" ) -0, nwi-1; swi-eye (1 +csc-nsc) ; end 
if exist (" swo" ) -0, nwo-1; swo-eye (1 + rsc-nsc) ; end 


if exist ("nwi" ) -0, di splay (" ERROR: nwi is not 
if exist ("nwo" ) -0, display (" ERROR: nwo is not 
(r, c) -size (swi) ; 

if r<>nwi+ma_tmp+me_tmp+lya+lye, ... 

display (" ERROR: number of rows in swi is not 
retf; end; 

if c < nwi, display (" ERROR: number of columns 
retf; end; 

(r, c) -size (swo) ; 
if c<>nwo+ka_tmp+ke_tmp, ... 
display (" ERROR: number of columns in swo not 
retf; end; 


defined in init .dat "); ret f; 
defined in init .dat "); ret f; 


consistant with other data");, 
in swi is less than nwi"); ... 


consistant with other data');. 



if r < nwo, display (' ERROR: number of rows in swo Is less than nwo' ) ; . .. 
retf; end; 

// Verify that the lamda provided is a row vector or a scalar and is 
// consistant with pea 

(rlambda, clambda ) -size (lambda) ; 

if rlambda > 1, di sp (' ERROR: lambda must be a row vector or a scalar!'), retf 
if clambda-1, lambda-lambda*ones (1, pea) , ... 
elseif clambda <> pea,... 

display (' ERROR: lambda must contain pea entries!'); retf; 

// The stopping conditions vector Is given a default value if it was not 

// provided (it will be saved in par.dat) 

if exist (' stop' ) -0, stop- (le-7; le-7; le-7; le2; le-2 ) ; 

save 'par.dat' stop 

// Inquire if the airframe or the engine submatrices are to be held fixed 
// storing the reply in variable aore 

// If not, put the subcontroller matrices in modified modal form. 

// Make ska and ske with orders nka and nke the variables to be used 

// from now on (Instead of s_ka and s_ke) 

nka-ns_ka; nke-ns_ke; 

inquire aore 'ENTER 1 TO FIX AIRFRAME, 2 TO FIX ENGINE or 0 FOR NEITHER: '; 
if aore-0, atot-nka +nke; ska-modi (s_ka, ns_ka) ; ske-modl (s_ke, nske) ; . . . 
elseif aore-1, atot-nke; ska»s_ka; ske-modl (s_ke, ns_ke) ;.. . 
elseif aore-2, atot-nka; ska-modi (s_ka, ns_ka) ; ske-s_ke;... 
else DISPLAY ('YOU MUST ENTER 0, 1 or 2!'); retf;... 
end, 

// Inquire if any of the D submatrices are to be fixed. DAA, DAYA, 

// DEE, and DEYE are determined by the corresponding submatrices of the 
// centralized controller. Under option 1 DEEA and DEAA are determined 
// by optimization. Under option 2 DEEA is sot to 0 and DEAA is determined 

// optimization. Under option 3 DEEA-0 and the remaining Ds are determined 

// by optimization. Under option 0 no Ds are fixed. 

DISPLAY (' ENTER 1 TO FIX ALL Ds EXCEPT DEAA t DEEA, 2 TO INCLUDE DEEA,'); 
inquire flxd ' 3 FOR ONLY DEEA, or 0 for NONE:'; 

if (fixd-1) Mfixd-2) Mflxd-3) *fixd <> 0,... 

DISPLAY ('YOU MUST ENTER 0, 1, 2 or 3!'); retf;end, 

// Get the row and column size of the plant matrix, input weighting matrix, 

// and the airframe and engine controller matrices. 

(rsp, csp) -size (sp) ; 

(rswi,cswi) -size (swi) ; 

(rska, cska) -size (ska) ; 

( rske, cskej -size (ske) ; 

// If the measurements from plant to controllers la,le, are absent, 

// create a dummy feedback of one (la - 1 and/or le - 1). 

// lya and lye hold the actual number of ya and ye measurements. 

// Increase the sizes of the system matrices to accomodate these 
// feedback variables. In the case of the dummy feedbacks being created 
// (and not input via init.dat) the system matrices are adjusted with 
// zero columns. 

If exist (' la' ) -0, ... 
la-1; lya-0; ... 
ska-(ska, 0*ones (rska, 1 ) J ; . . . 
if exist (' le' ) -0, .... 
le-1; lye-0; .. . 
sc-(sc, 0*ones (rsc, 2) ) ; . . . 
sp-(sp; 0*ones (2, csp) ) ; , . . 
swi-(swi; 0* ones (2, cswi )];... 
ske-(ske, 0*ones (rske, 1) );.. . 
else lye-le;... 
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sc- (sc ( : , (1 :csc-le) ) , 0*ones (rsc, 1) , sc (: , (csc-le+1 :csc) ) ) ; . . . 
sp-(sp( (1 :rsp-le) , :); O'ones |l,csp) ; sp ( ( rsp-le+1 : rsp) , :));... 
swi- (swi ((l:rswi-le|, :) ; 0*ones (1 , cswi) ; swi ( (rswi-le+1 :rswi) , :) );... 
end; . • . 
else ... 

lya-la; ... 
if exist (' le' ) -0, . . . 
le-1; lye-0; ... 
sc-(sc, 0*ones (rsc, 1) ) ; . . . 
sp-(sp; 0 # ones (1, csp) J ; ... 
ske-[ske, 0*ones (rske, 1) ) ; . .. 
swi-(swi; 0*ones (1, cswi) ) ; . . . 
else. . . 

lye-le; . . . 
end, . . . 
end, 

// Decompose the airframe subcontroller in modified modal form 

// ska- [aa, ba;ca, da ) . Use da to figure out the total Inputs matot 

// and outputs katot for ska. Subtract the la input from matot 

// to get ma, the actual number of inputs to ska. Subtract the number 

// of pea airframe to engine commands from the katot to get ka the number 

// of outputs going to the plant 

(aa ba ca da) -spl it (ska, nka) ; 

(katot matot ) -size (da) ; 

ma-matot-la; 

ka-katot-pea; 

// Decompose the engine subcontroller in modified modal form 

// ske-(ae, be;ce, de) . Use de to find ke the number of outputs 

// from the engine to the plant. Also find metot, the total number of 

// inputs, and subtratk pea, the airframe to plant commands, and le, the 

// engine feedback to get me, the actual number of external inputs 

// to the engine subcontroller 

(ae be ce de J -spl it (ske, nke) ; 

(ke metot ] -si ze (de) ; 
me-metot-pea-le; 

// If any of the entries of the first column of ba or be are <- zero (or 
// almost <- 0) then set them to le-9 (almost zero but positive). 

// This is required by the modified modal form of the system matrix, 
for i-1 :nka, if abs (ba ( i, 1) ) <le-9, ba (i , 1 ) -le-9; end, end, bal-ba(:,l); 
for i-1 :nke, if abs (be (i, 1) ) <le-9, be (i, 1) -le-9; end, end, bel-be(:,l); 

// Decompose the integrated plant sp- (ap, bp; cp, dp) 

(ap, bp, cp, dp) -split (sp, np) ; 

(kptot ntmp) -si ze (cp) ; 
kplnt-kptot-pea; 

(ntmp mptot )-size (bp) ; 
kpmpea-kptot-pea+1; 

// Let sk-sc be the centralized controller 
sk-sc; nk-nsc; 

// Decompose the centralized controller sk- (ak, bk; ck, dk | 

(ak,bk,ck,dk) -split (sk,nk); 

// If the user has not entered the (optional) matrix Nperf for 
// normalizing the performance cost, it is set to a matrix of ones 
// This matrix should be a frq(3,l)xl matrix, 
if exist ('Nperf ') -1, ... 

Nperf - max (Nperf , le-13) ; ( r, c) -size (Nperf );.. . 

if rofrq (3,1), ... 



display (' ERROR: The number of performance weights Is different from the ... 
number of frequency points ' ) ; retf ?end; . . . 
else, Nperf-ones (f rq (3, 1 ) , 1) ; . . . 
end; 

// If the user has not entered the (optional) matrix Ntrack for 
// normalizing the tracking cost, it is set to a matrix of ones. 

// This matrix should be a frq(3,l)xpea matrix, 
if exist (' Ntrack' ) -1, ... 

Ntrack - max (Ntrack, le-13) ; (r, c) -si ze (Ntrack) ;.. . 

If rofrq (3,1) , . . . 

display (' ERROR: The number of tracking weights disagrees with the ... 
number of frequency points ' ) ; ret f ;end; . . . 

1 f copea, ... 

display (' ERROR: The number of tracking weights disagrees with the ... 
number of interface variables ' ) ? retf ;end; . . . 
else, Ntrack-ones (f rq (3, 1 ) , pea) ; . . . 
end; 

// Define the weights for Simpsons rule integration. They depend on the 
// number of observation points specified in frq(3) and are stored in 
// the variable WEIGHT, 
weight- (1 ) ; 

for i-1: ( (frq (3,1) -3) /2) , weight- (weight , 4, 2 1 ; end 
weight- (weight, 4,1); 

// Define the set of frequency points at which the cost function and gradient 

// are evaluated. Store as the variable OMEGA.. 

omega -frq (1)* (( (frq(2)/frq(l))** (1/ (f rq (3) -1) ) ) * MO : f rq (3) -1 )) ? 

// Decompose the plant controller, splitting up the matrix bp relative to 
// its inputs and the matrix cp relative to its outputs. For example, 

// cap would be the plant to airframe piece, the first ma rows of the 
// plant matrix because these rows form the input to the airframe. 

(mp ntmp) -size (cp) ; 
kp- ka4ke; kapl- ka + 1; 

cap- cp([l:ma), :) ; ceap- cp ( ( (ma+1) : (ma+pea) ) , : ) ; 
cep- cp ( (ma+pea+1) : (ma+pea+me) , : ) ; 
cyap-cp(( (ma+me+pea+1) : (ma+me+pea+la) ), :) ? 
cyep-cp ( ( (ma+me+pea+la+ 1) : (ma+me+pea+la+le) ], :); 
bpa- bp(:, (l:ka)); bpe- bp(:, ( (ka+1) :kp) ) ; 

// The centralized controller is decomposed in a similar way as the plant 

// controller 

(kc ktmp) -size (ck ) ; 

bak- bk ( : , 1 :ma) ; bek- bk (: , ( (ma+1) : (me+ma) ) ) ; 

byak- bk (:, ( (ma+me+1) : (ma+me+la) ) ) ; byek-bk (:, ( (ma+me+la+1) : (ma+me+ la+le) )) ; 

cak- ck (1 :ka, :) ; cek- ck (kapl : kc, : ) ; 

dak- dk((ljka), (lsma)); daek-dk ((l:ka) , (ma+l:ma+me) ) ; 

dayek- dk ((l:ka), (ma+me+la+1 :ma+me+ la+le) ) ; 

deak-dk ( ( ka + 1 :ka+ke ) , (lima)) ;deyak-dk l (ka+1 :ka+ke) , (ma+me + 1 -.ma + me+la)) ? 

dek-dk ( (ka+l:ka+ke) , (ma+1 :ma+me) ) ; 

dyak-dk ((l:ka), (ma+me+1 :ma+me+ la )) ; 

dyek-dk ( ( ka+1 :ka+ke ) , (ma+la+me+1) : (ma+la+me+le) ); 

// HERE Construct the state space representation of Tcent from input: zac to 
// output: zea using the system with the global controller. This 

// is used in evaluating the zea tracking cost, but does not alter 
// with the choice of parameters so it is constructed once only 
// in this routine and then stored. 

ag- (ak, (-bak*cap-bek*cep+byak*cyap+byek*cyep) ; . . . 

(bpa*cak+ bpe*cek), (ap- (bpa*dak*cap) - (bpe*dek*cep) + . . . 

bpa*dyak *cyap+bpe*dyek *cyep-bpa*daek*cep+bpa *dayek *cyep. . . 

-bpe*deak*cap+bpe*deyak*cyap) ) ; 




.mtx 

bg- (bak; bpa* dak); 
eg- (0*ones (pea, nk) ceap); 
dg- 0*ones (pea, ma) ; 

// Let sg be the constructed state space representation of Tcent 
sg* [ag, bg; eg, dg) ; 
ng- nk+np; 

// If Ds are fixed, their initial values are determined here 

if aoreol, if (flxd-1) Mf lxd-2) -0,da-[dak,dyak;da ((ka+lska+pea) , :)| ;end;end; 

if aore <> 2, if (f ixd-2) * (f lxd-3) -0, de ( : , (1 :pea ) ) -0‘ones (Ice, pea) ; end;... 

if (f lxd-1) * lflxd-2) -0, de-(de(:, [ 1 :pea| ) , dek, dyekl ; end;... 

end; 

daa-da ( ( 1 : k a ) , (1 :ma ) ) ; daya-da ( (1 :ka) , (ma+1 :ma+la ) ) ; 
dee-de ( : , [pea+1 :pea+me) ) ; deye=de(:, (pea+me+1 :metot ) ) ; 
ska- ( aa , ba; ca, da); ske-(ae,be; ce, del; 



// Store the current (initial) values of ska and ske in skaO and skeO 
// to be used if the subcontrol ler (s) are fixed. 
skaO-ska; skeO-ske; 

// The stability upperbound stabil for the A submatrices should 
// be read in from init.dat and later stored in const.dat 
// If it is not defined or it is defined too high, it is given a 
// default value 

if exist (' stabil') -0 then stabil— le-9; 
if stabi l>-le-9, stabil — le-9; 

// Save the constants in const.dat 

save 'const.dat' la le lya lye ka ma ke me pea nka nke sp sk np nk frq ... 
omega daa dee daya deye atot sg ng swi nwi swo nwo weight stabil ... 
bal bel skaO skeO aore Nperf Ntrack fixd; 

// if there was no error in the data checking 
p_i-longcol (ska, ske) ; 

[ row, col ) -size (p_l) ; 
hi-eye (row) ; 
fx-0*ones (atot, 1) ; 

for 1-1 :atot, .... 

if p_i (I) >stabi 1, ... 
p_l (I)-stabil; ... 
fx(l) -1; ... 
hi (I, I) -0; . . . 
end; end; 

// Initialize costs and gradients 
(gradi,gradO,gradl, fh, jhO, jhl)-cost (p_i, lambda) ; 

// Save intermediate results 

save 'inter.dat' p_l gradi gradO gradl fh jhO j hi hi lambda fx; 

// Save starting data for restarts 

save 'start.dat' aore p_l gradi gradO gradl fh jhO jhl lambda fx; 

// call the PARTITION routine returning the final optimized vector 
p_opt -PARTITION (stop) ; 

// Generate the optimized ska and ske as the final output of the program 
( ska_opt , ske_opt) -MAT (p_opt) ; 
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load 'inter.dat' p_l fh jhO Jhl gradl gradO gradl hi fx lambda; 

// If a dummy feedback for ya or ye was added, then remove It 
If lya-0, (r, c) -si ze (ska_opt ) ; ska_opt-ska_opt (: , 1 :c-l ) ; end; 

If lye-0, (r,c) -si ze (ske_opt) ; ske_opt-ske_opt ( : , 1 : c-1 ) ; end; 

// save final data 

save ' inter. dat' p_i fh jhO Jhl gradl gradO gradl hi fx lambda ska_opt ske_opt; 


ret f 


// The UDF start. mtx is the main routine for the controller partitioning 
// program. It expects the input to the program to be In the datafile 
// 'init.dat', and requires certain data to be there. 

// 'init.dat' contains the global controller sc, the initial partitioned 
// controllers s_ka and s_ke, and the plant sp along with their orders nsc, 

// nska, nske, and np resp. It also contains the number of intermediate 
// variables pea and the range of frequencies for calculating norms frq. The 
// optional weighting matrices swl and swo (for Input and/or output weighting) 

// may be included (default - Id in both cases) as may la and/or le, the 
// numbers of neg. feedback variables (if these are absent, then the code 
// inserts dummy values of 1 and loads appropriate zeros In the data matrices.) 
// stabil - upper bound for eigenvalues of Aa and Ae to assure stability is 
// also contained in init.dat (if missing, defaults to -le-9) . 

// 

// The code gives the option of fixing values of either subcontroller If 
// optimization over the parameters in only one subcontroller is desired. 

// 

// The option is also given to fix certain of the D submatrices. 

// 

// The subcontrollers are put into a modified modal form where the A matrices 
// consist of 2x2 companion blocks and the first column of the B matrix is 
^ // fixed to its initial value (or le-9 if this initial value is zero) . These 
// are then split to identify the state space representations from inputs to 
// outputs using intermediate variables. 

// Here the modifications are made to account for the absence of la and/or le. 
// 

// The logarithmic frequency range over which the cost function is evaluated is 
// created as Omega. The vector of Weights used in Simpson's Rule Integration 
// is also calculated. 

// 

// The state-space representation for the nominal "command-t racking" 

// transfer matrix Tcent:zac to zea is computed, sg. Its row norms are stored 
// to be used as normalizations in the fTrack part of the cost function. 

// 

// The file 'const.dat' contains constants — - the numbers of subcontroller 
// inputs, ma and me; outputs, ka and ke; and intermediate variables, pea; 

// the numbers of negative feedbacks from the plant (lya/lye relects the actual 
// numbers, i.e. 0 If la resp le were absent In 'init.dat'. In this case la/le 
// would contain the dummy value 1) . Further constants are the state space 
// matrices for the plant, sp; the centralized controller, sk; the 
// nominal tracking transfer matrix, sg; the frequency weighting matrices, 

// swl and swo; and their orders np, nk, ng, and nw along with the orders 

// of the subcontrollers nka and nke; plus the frequency range, frq; the 
// vector of frequencies, omega; the weights for Simpson's Rule, weight; and 
// the vector of row norms for the nominal tracking transfer matrix. Other 
// constants are the values of the D subcontroller submatrices which are fixed 
// equal to the corresponding submatrices of the state space representation 
// for the centralized controller; the first columns of the subcontroller B 
// matrices which are fixed by the modi canonical form; and finally, the 
// values of skaO and skeO which are used In case the appropriate subcontroller 
// is fixed, atot Indicates the sum of orders of the subcontroller state-space 
// matrices which are not fixed. 

// 

// The initial system matrices for the subcontrollers skaO and skeO are put 
// into the form of a parameter vector p_l, an identity matrix Is created as 


// the Initial approximation for the Inverse Hessian used in the optimization 
// process and a vector, fx, indicating which of the entries in the state-space 
// A matrices are at the stability upper bound, pi and lambda are passed to 
// the COST function which returns the initial values of: 

// fh - the total cost 
// jho - the performance cost 
// Jhl - the tracking cost 

// gradl - gradient of the total cost 

// gradO - gradient of the performance cost 

// gradi - gradient of the tracking cost 

// 

// These values are stored in 'lnter.dat' the file for intermediate results. 

// 'start.dat' contains a record of the starting configuration. 

// 

// The parameters are passed to PARTITION which carries out the optimization. 

// This procedure reads the data It needs from the file 'inter.dat' and 
// returns the optimized paramter vector p_opt. 

// This final vector is split into the optimized airframe and engine 
// controllers ska_opt and ske_opt, which are the final output of 
// the program. 

// Created by Phil Schmidt and students Nader Kamranl and Brian Holawecky at 
// The University of Akron with the support of NASA Lewis Research Center 
// under grant NAG-3-1146. 
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//<zer> - z (nrow, ncol) 

// Construct a matrix of zeros of size nrow x ncol. This is done so that 
// less space is taken within the main program 
// Input: 

// nrow - number of rows 

// ncol - number of columns 

// output: 

// zer - the zero matrix of the desired size (nrow, ncol) 


zer-0*ones (nrow, ncol) ; 


ret f 

// Created by Phil Schmidt at The University of Akron with the support 
// of NASA Lewis Research Center under grant NAG-3-1146. 
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zero.mtx 



Appendix IV 
Data for Example 


65 


The data file IN1T.DAT corresponding to the example contains the following: 

PEA = 1 — number of interface variables 

FRQ = [0.1; 100; 41] — describes the frequency range of 41 points between 0.1 and 100 

NP =13 — order of the integrated plant 

NgC =13 — order of the centralized controller 

NS-KA =10 — order of the airframe subcontroller 

NS-KE = 7 — order of the engine subcontroller 

NWI = 26 — order of the input weighting transfer matrix 

NTRACK = 


6.1823D — 01 6.3153D - 01 

6.4195D — 01 

6.5131D — 01 

6.6139D — 01 

6.7411D — 01 • • 

6.9179D - 01 

7.1732D — 01 

7.5413D - 01 

8.0582D - 01 

8.7492D — 01 • • • 

9.6037D — 01 

1.0534D + 00 

1.1326D + 00 

1.1646D + 00 

1.1194D -t- 00 • • 

9.9744D — 01 

8.3547D - 01 

6.8055D - 01 

5.6432D — 01 

4.9513D — 01 • • 

4.6326D — 01 

4.5302D - 01 

4.5255D - 01 

4.5595D - 01 

4.6188D — 01 • • 

4.7282D — 01 

4.9653D - 01 

5.5002D - 01 

6.6494D - 01 

8.9018D — 01 • • 

1.2933D + 00 

1.9711D + 00 

3.0733D + 00 

4.8371D + 00 

7.6375D + 00 • • 
« 

1.2067D + 01 

1.9074D + 01 

3.0171D + 01 

4.7762D + 01 

7.5645D + 01 


41 x 1 vector of tracking weights (The size of NTRACK is FRQ(3)x PEA.) 

SP (listed below) — state-space matrix for the integrated plant transfer matrix 

SC (listed below) — state-space matrix for the centralized controller transfer matrix 

S_KA (listed below) — state-space matrix for the initial airframe subcontroller transfer matrix 

S-KE (listed below) — state-space matrix for the initial engine subcontroller transfer matrix 

SWI (described below) — state-space matrix for the input weighting transfer matrix transfer matrix 


The data matrices and initial partitioning matrices for the controller partitioning example are listed below. In all cases 
the given matrices A, B, C and D correspond to the state-space representation of the given system or subsystem 


^ = Ax + By 
d t 

u = Cx + Dy. 


A B' 
C D ' 


The state-space matrix for this system is 5 where 5 - 
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where 


The integrated airframe propulsion system with integrator augmentation is represented by SP - 


AP 

CP 


BP 

DP 


-4.40E' 2 

3.60E- 

-2 

-3.85E+ 1 


-3.18E+ 1 

1.40E' 2 

3.14 E~ 4 

2.59E -4 3.81E -2 


-2.27 E -1 

— 4.46E" 

-1 

1.94E+ 2 


-4.59E+ 0 

5.19E -4 - 

-1.57E -5 - 

2.10E -6 1.82 E~ A 


— 3.09J57- 3 

1.51E’ 

-2 

— 1.94E -1 


-4.81E- 

-4 

2.56E- 5 

9.46 E- 7 

3.74 E" 7 3.66E -5 


0 

0 


1.00E+ 0 


0 


0 

0 

0 

0 


1.42E -1 

— 9.89E" 

-1 
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The augmented plant state variables are 


x 


W,W,q 


,0,h,Ar2P,AT25,P6, 


T41B, 6tv> W F, >178, .A8] 


where 

u = aircraft body axis forward velocity, ft/s 
u; = aircraft body axis vertical velocity, ft/s 
q = aircraft pitch rate, rad/s 
6 = pitch attitude, deg 
h = altitude, ft 

A r 2F = engine fan speed, rpm 

N2S = core compressor speed, rpm 

P6 = engine mixing plane pressure, psia 

TA1B = engine high pressure turbine blade temperature, 0 R 
<5 T v = normalized thrust vectoring angle, deg/ 10 
WF = normalized engine main burner fuel flow rate, lb m/hr/500 
A78 = normalized thrust reverser port area, in 2 /50 
AS = normalized main nozzle throat area, in 2 / 100. 


The augmented plant outputs are listed in the order 


’ 2a ’ 

with z a — 

v m 


N2 ‘ 

2 e 

Qv 

, 2 C — 

EPR_ 







and z M = FEX where 


V = normalized aircraft speed, ft/s/20 

q v = normalized pitch variable, (g(deg/s)-f O.10(deg))/3 

N2 = normalized engine fan speed, % of maximum allowable rpm at operating condition/5 


EPR = normalized engine pressure ratio, ratio/0.3 

The augmented plant inputs are the rates of change on normalized control variables as described in the text. 
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The centralized controller has state-space representation SC = 


AC BC 
CC DC 


where 
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[ A.KA B_KAl , 

The initial partitioning consisted of the airframe controller with representation S-kA - I D_KAJ ° order 1U 

and an engine controller with representation S.KE = of order 7. The initial partitioning was obtained by 

the procedure of reference [ 2 ]. This partitioning is not in the minimal parameter form. The submatrices in the initial 
partitioning are: 
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Recall that the input weighting transfer matrix is VFi(s) = G(s)(J+ A(s)G(s)) 1 where G(s) and A(s) are the integrated 
plant and centralized controller transfer matrices respectively. The state-space representation for this transfer matrix is too 
large to list. Instead, the user can easily construct it from the following MATRJX x command applied to the state space 
matrices for the integrated plant and the centralized controller 

[SWI,NWI]=FEEDBACK(SC,NSC,SP,NP). 

Notice that this command also produces the correct value for NWI, the order of the weighting transfer matrix. 


The parameter optimization algorithm for controller partitioning was applied to the problem with the initial partitioning 
given by S_KA and S_KE as listed above and with input weighting SWI as described above. The controllers obtained from 
this process had state-space representations SKA.OPT for the airframe and SKE.OPT for the engine where, as before 


A-A.OPT = 


0 

-4.43E" 2 

0 

0 

0 

0 

0 

0 

0 

0 


1.00E+ 0 

-3.15E" 1 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

-7.28E+ 1 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 



0 

0 

0 

0 



1.00E + ° 

0 

0 

0 



1.12E+ 1 

0 

0 

0 



0 

0 

1.00E +o 

0 



0 

-1.72E+ 3 

-5.9IE+ 1 

0 



0 

0 

0 

0 



0 

0 

0 

-1.00E- 

■s 


0 

0 

0 

0 



0 

0 

0 

0 





0 

0 

0 



0 

0 

0 



0 

0 

0 



0 

0 

0 



0 

0 

0 



0 

0 

0 



1.00E +o 

0 

0 



-3.41E+ 0 

0 

0 



0 

0 

1.00E +o 


0 -9.96E" 3 — 1.99E -1 . 


71 


B-A.OPT = 
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The optimized engine controller has state-space representation 
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Residualization of high frequency modes was applied to the optimized engine controller to reduce it to one with order 4 
(not shown here because it is easily obtained). Balanced model reduction was applied to the optimized airframe controller 
to reduce it to one of order 6. The optimization procedure was applied to this sixth order subcontroller with the engine 
controller fixed at the one of fourth order. The resulting reduced order optimized airframe subcontroller is 
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